Computer  Modeling  of  Ceramic  Boride  Composites 


Dr.  Valeriy  V.  Kartuzov 

SCIENCE  AND  TECHNOLOGY  CENTER  IN  UKRAINE 
METALISTIV  7A,  KYIV,  UKRAINE 

INSTITUTE  FOR  PROBLEMS  OF  MATERIALS  SCIENCE 
3  KRZHYZHANOVSKY  STR.,  KYIV,  03680  UKRAINE 


EOARD  #STCU  P-510/STCU  11-8003 


Report  Date:  November  2014 
Final  Report  from  1  November  201 1  to  31  Oetober  2014 


Distribution  Statement  A:  Approved  for  public  release  distribution  is  unlimited. 


Air  Force  Research  Laboratory 
Air  Force  Office  of  Scientific  Research 
European  Office  of  Aerospace  Research  and  Development 
Unit  4515,  APO  AE  09421-4515 


REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  Information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  the  burden,  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply 
with  a  collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  fDD-MM-yyyy;  2.  REPORT  TYPE  3.  DATES  COVERED  (From  -  Toj 

01  November  2014  Final  Report  1  November  2011  -31  October  2014 


4.  TITLE  AND  SUBTITLE  5a.  CONTRACT  NUMBER 


Computer  Modeling  of  Ceramic  Boride  Composites 


STCU  P-510 


5b.  GRANT  NUMBER 


STCU  11-8003 


5c.  PROGRAM  ELEMENT  NUMBER 


6.  AUTHOR(S) 


61102F 

5d.  PROJECT  NUMBER 


Dr.  Valeriy  V.  Kartuzov 


5d.  TASK  NUMBER 


5e.  WORK  UNIT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

SCIENCE  AND  TECHNOLOGY  CENTER  IN  UKRAINE 
METALISTIV  7A,  KYIV,  UKRAINE 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


INSTITUTE  FOR  PROBLEMS  OF  MATERIALS  SCIENCE 
3  KRZHYZHANOVSKY  STR.,  KYIV,  03680  UKRAINE 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

EOARD 
Unit  4515 

APO  AE  09421-4515 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 
AFRL/AFOSR/lOE  (EOARD) 


11.  SPONSOR/MONITOR’S  REPORT  NUMBER(S) 
AFRL-AFOSR-UK-TR-201 5-001 6 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


13.  SUPPLEMENTARY  NOTES 


14.  ABSTRACT 

This  effort  was  to  investigate  regularities  and  mechanisms  of  kinetics  growth  of  boride-boride  composites  of  eutectic  co¬ 
crystallization  by  means  of  computer  modeling.  In  terms  of  project's  scope  of  work  a  lot  of  new  important  results  were  obtained. 
Calculations  showed  that  coefficients  of  thermal  expansion  (CTE)  for  LaBs  and  systems  LaBs  -  MeB2  (Me  -  Zr,  Ti)  have  an 
jump  in  temperature  interval  from  750  K  to  1 000  K.  At  elevation  of  temperature  the  energy  of  contact  surface  in  these  systems 
decreases  however  the  system  preserves  its  mechanical  properties  up  to  melting  temperatures.  Ideal  work  of  adhesion  of  the 
interface  [001]-LaB6//[001]-TixZri.xB2,  (1 10)-LaB6//(110)-TixZri.xB2  increases  if  titanium  concentration  elevates.  Change  of  a 
direction  of  fibers  growth  at  directed  solidification  of  eutectics  LaB6-TixZri-xB2  is  associated  with  a  big  loss  in  energy  at 
formation  of  semi-coherent  interface  as  titanium  concentration  increases.  A  conception  of  internal  energy  sources  allocated  on 
interfaces  and  in  phase  volume  of  composites  and  is  the  bond  energy  among  composite's  components  in  a  solid  state:  surface 
energy  of  interfaces,  spot  defects,  internal  stresses  governed  by  components  CTE,  etc.  is  proposed.  A  new  model  of  interface 
evolution  at  heating  where  a  surface  boundary  between  phases  turns  into  a  boundary  with  a  thickness  (i.e.  volumetric)  and 
consists  of  two  solid  parts  with  a  liquid  between  is  presented.  Computer  software  to  obtain  multi-fractal  characteristics  of  2D 
images  obtained  in  electron  microscope  is  developed.  A  package  of  computer  programs  for  2D  imitative  modeling  of  formation 
of  eutectic  structures  is  developed  on  the  base  of  cellular  automates.  Moreover  the  analytical  review  of  both  a  retrospective  and 
current  Russian  and  Ukrainian  papers  in  the  field  of  new  materials  for  aviation  is  carried  out.  The  results  obtained  under  this 
project  are  published  in  21  papers,  one  monograph  is  also  submitted  to  print;  17  papers  are  presented  at  conferences. 


15.  SUBJECT  TERMS 


EOARD,  Materials,  mictrostructural  characterization,  high  temperature 


17.  LIMITATION  OF  18,  NUMBER 
ABSTRACT  OF  PAGES 


16.  SECURITY  CLASSIFICATION  OF: 

a.  REPORT 

UNCLAS 

b.  ABSTRACT 

UNCLAS 

c.  THIS  PAGE 

UNCLAS 

19a.  NAME  OF  RESPONSIBLE  PERSON 

Matthew  Snyder 

19b.  TELEPHONE  NUMBER  (include  area  code) 

+44  (0)1895  616420 


Standard  Form  298  (Rev.  8/98) 

Prescribed  by  ANSI  Std.  Z39-18 


Project  Title:  COMPUTER  MODELING  OF  CERAMIC  BORIDE 

COMPOSITES 


Final  Report 

Reporting  Period:  01/11/2011  -  01/11/2014 


Principal  Investigator:  Kartuzov  Valeriy  Vasilievich  (PhD) 


Institution:  Frantsevich  Institute 
for  Problems  of  Materials  Science 
of  National  Academy 
of  Sciences,  Ukraine 


FOARD  Project  Number:  118003 
STCU  Project  Number:  P5  10 


Distribution  A:  Approved  for  pubiic  reiease;  distribution  is  uniimited. 


ABSTRACT 

This  effort  was  to  investigate  regularities  and  meehanisms  of  kineties  growth 
of  boride-boride  eomposites  of  euteetie  eo-erystallization  by  means  of  eomputer 
modeling.  In  terms  of  projeet's  seope  of  work  a  lot  of  new  important  results  were 
obtained.  Caleulations  showed  that  eoeffieients  of  thermal  expansion  (CTE)  for 
LaBe  and  systems  LaB6  -  MeB2  (Me  -  Zr,  Ti)  have  an  jump  in  temperature  interval 
from  750  /if  to  1000  K.  At  elevation  of  temperature  the  energy  of  eontaet  surfaee  in 
these  systems  deereases  however  the  system  preserves  its  meehanieal  properties  up 
to  melting  temperatures.  Ideal  work  of  adhesion  of  the  interfaee  [001]-LaB6//[001]- 
TixZri.xB2,  (110)-LaB6//(110)-TixZri.xB2  inereases  if  titanium  eoneentration 
elevates.  Change  of  a  direetion  of  fibers  growth  at  direeted  solidifieation  of 
euteeties  LaB6-Ti^Zri_xB2  is  assoeiated  with  a  big  loss  in  energy  at  formation  of 
semi-eoherent  interfaee  as  titanium  eoneentration  inereases.  A  eoneeption  of 
internal  energy  sourees  alloeated  on  interfaees  and  in  phase  volume  of  eomposites 
and  is  the  bond  energy  among  eomposite's  eomponents  in  a  solid  state:  surfaee 
energy  of  interfaees,  spot  defeets,  internal  stresses  governed  by  eomponents  CTE, 
ete.  is  proposed.  A  new  model  of  interfaee  evolution  at  heating  where  a  surfaee 
boundary  between  phases  turns  into  a  boundary  with  a  thiekness  (i.e.  volumetrie) 
and  eonsists  of  two  solid  parts  with  a  liquid  between  is  presented.  Computer 
software  to  obtain  multi-fraetal  eharaeteristies  of  2D  images  obtained  in  eleetron 
mieroseope  is  developed.  A  paekage  of  eomputer  programs  for  2D  imitative 
modeling  of  formation  of  euteetie  struetures  is  developed  on  the  base  of  eellular 
automates.  Moreover  the  analytieal  review  of  both  a  retrospeetive  and  eurrent 
Russian  and  Ukrainian  papers  in  the  field  of  new  materials  for  aviation  is  earried 
out.  The  results  obtained  under  this  projeet  are  published  in  21  papers,  one 
monograph  is  also  submitted  to  print;  17  papers  are  presented  at  eonferenees. 
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INTRODUCTION 


This  effort  was  to  investigate  regularities  and  meehanisms  of  kineties  growth  of 
boride-boride  eomposites  of  euteetie  eo-erystallization  by  means  of  eomputer 
modeling.  The  approach  used  is  based  on  development  of  mathematical  models 
which  reliably  describe  this  process  through  the  whole  material  structure  scales 
(macro,  meso  and  micro)  and  conduction  of  numerical  experiment  targeted  on 
optimization  of  structure  and  service  properties  of  the  above  said  composites.  In 
terms  of  project's  scope  of  work  in  a  triad  "Model  (B)  -  Computer  program  (C)  - 
Result  of  numerical  experiment  (A)"  many  new  significant  results  were  obtained. 
The  most  striking  and  interesting  are  the  following: 

A.  (Results) 

>  Calculations  on  the  base  of  the  method  of  apriori  pseudopotential  showed 
that  coefficients  of  thermal  expansion  (CTE)  for  LaBa  and  systems  LaBe  -  MeB2 
(Me  -  Zr,  Ti)  have  an  jump  in  temperature  interval  from  750  K  to  1000  K.  At 
elevation  of  temperature  the  energy  of  contact  surface  in  these  systems  decreases, 
however  the  system  preserves  its  mechanical  properties  up  to  melting  temperatures 
(eutectic  temperatures).  The  presence  of  local  maximum  of  strength  is  shown  for 
the  composite  in  the  eutectic  point. 

>  Ideal  work  of  adhesion  of  the  interface  [001]-LaB6//[001]-TixZri.xB2,  (110)- 
LaB6//(110)-TixZri.xB2  increases  if  titanium  concentration  elevates.  Change  of  a 
direction  of  fibers  growth  at  directed  solidification  of  eutectics  LaB6-TixZri_^B2  is 
associated  with  a  big  loss  in  energy  at  formation  of  semi-coherent  interface  as 
titanium  concentration  increases. 

>  Modeling  of  interface  phenomena  in  Cu-Si  predicts  the  existence  of  a  new 
type  of  equilibrium  states  (solid  phase-liquid-solid  phase)  at  temperature  above  the 
eutectic  one;  these  states  describe  a  complete  disintegration  with  a  liquid  interface. 
This  shows  the  role  of  interface  energy  as  possible  mechanism  for  contact  melting 
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phenomena  and  formation  of  diffusion  zone.  The  obtained  results  also  explain  a 
phenomenon  of  separation  phases  in  liquid  eutectic. 

B.  (Models) 

>  A  conception  of  internal  energy  sources  allocated  on  interfaces  and  in  phase 
volume  of  composites,  which  reflect  the  energy  of  binding  between  composite's 
components  in  a  solid  state:  surface  energy  of  interfaces,  point  defects,  internal 
stresses  due  to  difference  of  CTE  of  components,  etc.  is  proposed.  Density  of 
internal  energy  sources  (it  depends  on  composite  structure  in  a  solid  state)  may 
cause  self-heating  of  composite  and  fall  of  energy  required  for  melting  of  the 
whole  composite,  i.e.  reduction  of  total  melting  temperature  of  composite  as  it 
occurs  in  eutectic  composites.  This  temperature  fall  can  be  significant  and  it 
depends  on  a  ratio  between  concentrations  of  composite  components  as  well  as 
structure  of  interphase  boundaries. 

>  A  new  model  is  presented  of  interface  evolution  at  heating  where  a  surface 
boundary  between  phases  turns  into  a  boundary  with  a  thickness  (i.e.  volumetric) 
and  consists  of  two  solid  parts  with  a  liquid  between.  The  model  is  developed  in 
assumption  that  diffusion  does  not  influence  temperature  distribution,  which  is 
defined  by  a  solution  of  separate  problem,  but  temperature  affects  diffusion 
through  the  dependence  of  diffusion  coefficient  on  temperature. 

C.  (Programs) 

>  Computer  software  to  obtain  multi-fractal  characteristics  of  2D  images 
obtained  in  electron  microscope  is  developed. 

>  A  package  of  computer  programs  for  2D  imitative  modeling  of  formation  of 
eutectic  structures  is  developed  on  the  base  of  cellular  automates  method. 

Moreover  the  analytical  review  of  both  a  retrospective  and  current  Russian 
and  Ukrainian  papers  in  the  field  of  new  materials  for  aviation  is  carried  out. 
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The  results  obtained  under  this  project  are  published  in  21  papers,  one 
monograph  is  also  submitted  to  print;  17  papers  are  presented  at  international 
conferences;  two  scientists  involved  (S.  Ivanov  and  A.Khachatrian)  defended  their 
PhD;  students  partially  involved  in  this  project  defended  4  Master's  diploma  works. 
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I.  MACROSCALE 


1.  Continuum  model  of  formation  of  eutectic  compositions  by  directional 
solidification  method  (macro  level). 


A  continuum  model  is  proposed  of  eutectic  ceramic  composites  formed  from 
two-component  melt  during  its  crystallization  (fig.l)  under  conditions  of  controlled 
heat  removal,  characterized  by  the  temperature  gradient  and  the  velocity  of  the 
crystallization  front  [1  -  9]. 


Fig.  1.  TEM  data  from  the  transverse  section  of  a  LaB6-ZrB2  DSE  (g)  is  the  SAD 
pattern,  which  illustrates  that  [0001]  ZrB2  is  approximately  parallel  to  [0001]  LaB6 
and  (110)  LaB6  parallels  ( 1100)  ZrB2.  Kikuchi  lines  evident  in  (c)  show  that  there 
is  actually  a  small  misorientation  between  two  phases.  Electron  diffraction  patterns 
from  individual  fibers  (b,  e,  f)  indicate  that  all  of  the  ZrB2  fibres  are  locally 
oriented  in  the  same  direction  within  0.02  degrees. 
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Figure  2  provides  examples  of  investigated  sehemes  of  regular  two- 
eomponent  euteetie  eomposites,  where  A,  B  —  mutually  insoluble  in  solid  state 
components  of  the  composite;  A  —  matrix  (e.g.,  LaBa),  B  —  fibers  with  diameter  d 
or  plates  with  thickness  d  (eg,  TiB2,  ZrB2,  HfB2,  ...).  It  is  assumed  that  composite 
component  B  has  eutectic  concentration  c,  see  Figure  2. 


Fig.  2.  Scheme  of  eutectic  compositions  with  regular  structure:  {a)  -  fibrous;  {b)  - 
lamellar. 

Formulation  of  the  initial-boundary  problem  of  directional  solidification  with 
explicit  isolation  of  the  interphase  interface. 

Fig.  3  presents  a  diagram  of  temperature  distribution  T{x,t)  in  the  melt 
(0  <  V <  S{t))  and  regular  composite  (in  solid  phase,  S{t)  <x<L)  under  non¬ 
stationary  mode  of  its  formation.  Accepted  notation:  v  =  const  —  velocity  of 
pulling  of  the  composite  in  direction  of  v  axis,  t  —  time,  Tc  —  temperature  of 
crystallization  and  structure  formation  of  the  composite. 

It  is  assumed  that  the  insulation  of  the  lateral  surface  is  provided  as  a  part  of 
equipment  design  and  lateral  (in  the  direction  perpendicular  to  the  axis  x)  heat 
removal  is  absent.  Its  presence  is  undesirable  because  it  facilitates  formation  of 
heterogeneity  of  composite  in  the  transverse  direction  (along  the  front  of 
crystallization).  Therefore,  heat  flow  takes  place  in  the  direction  of  the  axis  v,  and 
depends  on  the  temperature  distribution  along  this  axis  (see  Fig.  3) 

The  problem  of  modeling  the  composite  formation  is  formulated  as  the 
following  modified  Stefan  problem: 
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Find  temperature  distribution  T[x,t): 


T  =  T(x,t) 


T^{^x,t),  ^<x<S{t) 
T^{^x,t),  S{t)<x<L 


(1) 


and  law  of  motion  of  front  of  erystallization  and  structure  formation  S(t). 


Fig.  3.  Scheme  of  temperature  distribution  T(x,t),  T^<T^,  t  >  0 . 

Fig.  4  provides  scheme  of  phase  diagram  of  the  eutectic  composite. 

T 

-'a 

T  =T  =T 


0  c  1 

(^)  (5) 

Fig.  4.  Schematic  phase  diagram  of  the  composite.  Tm  —  melting  temperature; 
Tc —  temperature  of  crystallization  (solidification);  Tg,  c  —  eutectic  temperature 
and  concentration  respectively. 

It  is  assumed  that  the  densities  of  the  liquid  and  solid  phases 
coincide:  pi  =  Ps  =  P-  From  this  and  from  the  continuity  of  flow  (p/V;  =  p^v^)  follows 
that  vi  =  Vs  =  V.  In  addition,  it  is  assumed  that  the  insulation  of  side  surface  of  the 
sample  is  provided  by  the  process  equipment  design  and  lateral  (in  the  direction 
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perpendicular  to  the  axis  x)  heat  flux  is  absent.  Its  presence  is  undesirable  because 
it  creates  non-uniformity  in  composite  in  the  transverse  direction  (along  the  front 
of  crystallization).  Because  of  this  heat  flow  takes  place  in  the  direction  of  the  axis 
X  and  depends  on  the  temperature  distribution  (1)  along  this  axis  (see.  Fig.  3), 
which  is  determined  by  the  equations: 


pc, 


pc. 


dT,_  dX 


dt 

dT^ 


dx^ 

dx^ 


pc,v 


pc.v 


dx’ 
dx 


0<x<S(f) 
S  (t)  <  X  <  L 


(2) 


Equations  (2)  are  complemented  by  the  boundary  conditions: 

A)  At  the  ends  of  interval  (0,L)  constant  temperature  is  maintained: 


7;(0;()  =  7;>7;,,  T,{L-,t)  =  T,<Z,  t>0. 


(3) 


B)  Conditions  at  the  front  of  crystallization  and  structure  formation  S(t): 

T,(S{t),t)  =  T,(S{t),t)  =  T,  t>0,  (4) 


dT, 


dx 


x=S-0 


dx 


+  {cs-Ci)pvT^=(5S{t),  S{t)<0,  (5) 


^=s+o 


0  =  0, +o^, 


(6) 


where  o  -  specific  (per  unit  volume)  energy  of  crystallization  of  composite 
components  and  formation  of  its  structure  (bonds  (contacts)  between  its 
components).  Specific  energy  o  has  two  components:  Oc  —  effective  latent  heat  of 
crystallization  (melting)  of  composite  per  unit  of  its  volume;  Ob  —  specific  energy 
required  for  the  formation  (breakup)  of  bonds  between  the  components  of  the 
composite. 
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Latent  heat  of  erystallization  per  unit  volume  of  the  eomposite  is  estimated  by 
the  formula 

(7) 

where  =  (l-c)^^  —  is  the  same  but  per  unit  mass;  i  —  latent  heat 

of  erystallization  (melting)  of  unit  mass  of  eorresponding  eomponents  A  and  B  of 
eomposite. 

For  energies  of  formation  (break  up)  of  boundaries  of  eomposites  {a)  and  {b), 
see  Figure  1,  related  to  unit  volume  of  the  eomposite  we  have,  respeetively, 
estimates 


dye  2yc 

a  a 


(8) 


wher  y  —  surfaee  energy  of  the  interfaee  between  eomponent  A  and  B. 

4c 

Composites  {a)  and  {b),  respeetively,  eontain  in  a  unit  volume  N  =  — 5-  of  parallel 

nd 


fibers  with  diameter  d  and  N  =  —  lamellae  of  thiekness  d  with  eoneentration  c. 

d 


From  relations  of  linear  thermoelastieity  we  obtain  the  following  assessment 
of  internal  stress  energy  per  unit  volume  of  euteetie  due  to  the  difference  of  the 
coefficients  of  thermal  expansion  of  the  composite  components 


G 


bl 


(9) 


where  E^={l-c)E^  +  cE^ ;  (3  =  £  +  (1  -  c)a,^E^  +  ca^E^^  AT;  E^,Eg  — 


moduli  of  elasticity  of  composite  solid  components;  oc^,  oc^  —  coefficients  of 
linear  thermal  expansion  of  composite  components;  AT  —  temperature  difference 
that  led  to  the  emergence  of  internal  stress  in  the  composite. 

Thus,  in  (6)  for  component  G^,  which  takes  into  account  bond  formation 
between  components  A  and  B,  we  have 


^b  =  ^bi 


+  G 


bl  ■ 


(10) 
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Conditions  (5),  (6)  differ  from  conditions  at  the  crystallization  front  of 
classic  Stefan  problem  because  they  take  into  account  not  only  crystallization  but 
also  structure  formation  of  composite,  component  o/,,  see  (10).  This  component 
differs  by  surface  energy  of  boundaries  between  composite  components  per  unit 
volume  of  solid  phase  o^i  and  by  energy  of  internal  stresses  (Jbi- 

Crystallization  temperature  or  eutectic  temperature  at  the  front  of 
crystallization  (see  Fig.  4)  is  given  by 

Cb9bC  +  C^^^{1-c) 


where  c^,Cj^  —  specific  (per  unit  mass)  heat  capacity  of  the  composite 

components.  Note  that  this  temperature  depends  on  the  structure  of  eutectic  and 
contact  interaction  between  its  components. 

Equations  (2)  are  complemented  also  by  initial  conditions 

T,(xfi)  =  T,(x)  =  T,,0<x<L-.  T,(Lfi)  =  T, 

T,(S,)  =  T,(S,)  =  T^,  S(r  =  0)  =  S„=L-n  ’ 


where  p  —  sufficiently  small  number,  because  at  the  initial  time  moment 
S  (O)  =  -oo ,  see.  Fig.  5. 

Technologically  important  is  steady  mode  of  pulling  the  sample  when  the 
sample  is  moved  with  a  constant  velocity  v.  The  conditions  of  front  stay  in  the 
selected  coordinate  system  and  passing  to  the  stationary  thermal  operation  mode 

( —  =  0)  are  as  follows: 
dt 

\imS(t)  =  S*  =  S*(v,d)  =  const  <L  (stationary  state  of  the  front),  (13) 
lim^  (t)  =  -v  (condition  of  the  front  stay) ,  (14) 


Distribution  A:  Approved  for  pubiic  reiease;  distribution  is  uniimited. 


=  0,  S*<x<L  (condition  of  stationarity).  (15) 


dT, 


dt 


=  0,  0<x<S*- 


dt 


Function  S(t)  is  decreasing  (left-hand  part  of  equation  (5)  —  negative),  and 
veloeity  5'(t)<0  —  inereasing  funetion.  The  front  stops  at  time  moment  t  =  t^ 

when  5  =  tp )  =  -V .  Beeause  of  this  equalities  (13),  (14)  take  the  form 

S{t  =  t^)  =  S{v,d,t^)  =  S*{v,d),  (16) 


S(t  =  t^)  =  S(v,d,t,)  =  -v . 


(17) 


Equation  (17)  is  the  equation  to  determine  the  time  moment  t  =  ty  when  the  front 
stops  and  steady  state  mode  begins. 

The  values  S*,  L  -  S*  should  signifieantly  exeeed  the  thiekness  of  the 
erystallization  front  A  □  .  Then  maero  formulation  of  the  problem  (when  front 


is  eonsidered  as  a  surfaee)  makes  sense. 

Equalities  (16),  (17),  in  essenee,  are  prerequisite  for  building  a  dependeney  d 


=  d(v).  This  relationship  exists  if  there  is  an  interval 


U<t^<t 


of  suitable  times 


t  =  t^,  that  is 


Eigure  5  shows  a  sehematie  solution  of  the  formulated  Stefan  problem.  The 
value  S*  indieates  the  position  of  the  front  of  erystallization  and  strueture  formation 
when  it  stopped,  that  is  when  S{t  =  t^)  =  -v .  The  eondition  of  rest  of  the  front  in 

the  ehosen  eoordinate  system:  =  tj-l-v  =  0.  The  initial  interval  of  funetional 

dependenee  S{t)  ean  be  approximated  by  the  dependenee  S{t)  =  L- ayft . 
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If  stationary  front  position  S(f^)  =  S*  is  known  from  the  solution  of  the  Stefan 
problem,  see  (16),  then  temperature  distribution 


T  =  T{x)^ 


T,(x),  ()<x<S* 
T^(x),  S*<x<L 


eorresponding  to  a  stationary  regime  of  erystallization  and  strueture  formation  of 
euteetie  composition  is  given  in  a  closed  form  by  formulas: 


T  -T 

T’  _ T’  I  ^0 


1-exp 


V 


5* 


y‘^1  j 


f 

f  \ 

f  W 

V 

V  . 

exp 

—  V 

-exp 

—S 

V 

J 

Jy 

,  a, 


\ 

pC/ 


,  0<v<5*,  (18) 


T'c-T'l 


1-exp 


''(L-V) 


f 

exp 

^x-S■) 

-exp 

1 

1 

* 

V 

(19) 


a  = — S  <x<L. 


For  temperature  gradients  at  the  front  x  =  S*  {v,d)  from  (18),  (19)  we  have 


dT, 

To-T,  V  1 

V  * 

dT^ 

> 

1 

1 

dx 

x=S*-0 

1-exp 

Vs- 

J 

- 

—  O 

) 

’  dx 

x=S*+0 

1-exp 

) 

From  this,  for  the  heat  flux  W  spent  on  this  front  for  crystallization  and  structure 
formation  of  eutectic  composite,  we  obtain 


W=x, _ - 

l-exp(v5yaj)«, 


exp 


-S 

\<^i  j 


—  Xc 


T.-T, 


l-exp(v(L-5*)/«,)«, 


-+(cj-c>vr,(20) 
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Fig.  5.  Schematic  solution  of  modified  Stefan  problem  and  initial  conditions, 
or 


c,{T„-T,) 

l-exp^vS'Ya,) 


exp 


v«/ 


(7; -71) 


l-exp(v(L-5'Y/a,) 


+  {c,-Ci)T^ 


pv.  (21) 


Conditions  (5)  at  the  front  x  =  S*(v,d)  take  form 

ty  =  ov  =  (o^  +  oJv 


(22) 


and  determine  the  implicit  functional  dependence  d(v)  of  structural  parameter  d 
from  velocity  of  directional  solidification  v.  For  example,  in  the  case  of  fiber 
composite  (a).  Fig.  2,  this  dependence  is  written  as 


l-exp(v<5'ya;) 


exp 


—S 

\^i  J 


c(7;-7l) 


l-exp(v(L-S'*)/«,) 


+{c,-c,)T=q+^+^ ,  (23) 


Jp  p 
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where  the  temperature  Tc  depends  on  the  fiber  diameter  d  and  the  size  S* {v,d)  of 

liquid  phase  depends  on  the  velocity  v  and  the  fiber  diameter  d.  For  lamellar 
composite  (b),  Fig.  2,  multiplier  factor  2  in  the  right-hand  part  of  (23)  should  be 
replaced  with  4. 


Calculations 

The  calculation  of  parameters  of  material  and  processes  was  carried  out  for 
LaB6-ZrB2  composite. 

Main  initial  data  are  given  in  the  Table  1. 

Table  1.  Properties  of  materials  of  the  composite. 


Material 

E, 

GPa 

a  (GTE) 

A 

kg/m^ 

Tm, 

°C 

kJ/ 

kg 

Csi 

J/ 

kg- grad 

Cl, 

J/ 

kg- grad 

A, 

w/ 

(m-grad) 

2/, 

W/ 

(m-grad) 

LaBe  [10] 

478.73 

13.24-10 

4720.0 

2210.0 

1798 

572 

1000’ 

47.7 

20’ 

ZrB2  [11] 

495.81 

28.0-10“^ 

6200.0 

(6090) 

3225.0 

930 

429 

953 

83 

58 

*)  approximate  estimate 


Thermal  properties  of  LaB6-ZrB2  composite  as  a  homogeneous  material,  such 
as  heat  capacity,  thermal  conductivity  and  specific  crystallization  (melting)  heat  for 
solid  and  liquid  phases  necessary  for  calculations  were  found  by  mixture  rule: 

Ci  =  (1  -  Co)c^LaB6  +  CffCsZvm\ 

Q  =  (1  -  Co)C/LaB6  +  Co-C/ZrB2; 

(1  -  Co)^sLaB6  +  Co-A,sZrB2', 

^1=  -  Co)Z/LaB6  +  Co-Z/ZrB2; 

^  =  (1  -  Co)^LaB6  +  Co-^ZrB2- 
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Table  2.  Calculated  properties  of  LaB6+ZrB2  composite 


Material 

q. 

Cl, 

c, 

h. 

kJ/ 

J/ 

J/ 

w/ 

W/ 

kg 

kg- grad 

kg- grad 

(mgrad) 

(m-grad) 

LaB6+ZrB2 

1616 

542 

990 

53 

26 

Other  initial  parameters  for  calculations: 

Co  =  0.16  (voL),  0.21  (mass);  y=  170.0  J/m^;  Tqq  =  20  °C;  Ty  =  0.9  Tm  (approximate 
estimate);  L  =  5  mm;  Tq=Tm+  100  °C;  Tl=Tm-  200  °C;  v  =  4.2-10”^  -  LS-IO"^ 
m/s. 


Fiber  diameter 


4yco 


IT 


CP(^c-^l) 


+ 


W-pq-w 

Cip(To-Tc) 


l-exp(v(L-5)/aJ  l-exp(v5/ aj 


-exp(v5/a;)  +  (c5  -c,)pTy., 


AC  (2;  =  ^  ,  (2^  =  ^  . 

PC  PC 


Table 


3.  Calculations  of  fiber  diameter  by  the  formula  (23)  {T^  =  2210  °C). 


V,  m/s 

5,  mm 

w,  J/m"' 

d,  |im 

110“^ 

1.25  (L=  5  mm) 

2.54-10’ 

1.1-10“^ 

2-10^^ 

2.5 

2.54-10’ 

6.7-10^^ 

5-10“^ 

1.25 

2.54-10’ 

1.28- 10“’ 

I-IO" 

1.25 

2.54-10’ 

-4.7-10“’ 

1-10  " 

3.75 

2.54-10’ 

1.6-10  ’ 

5-10“^ 

3.75 

2.54-10’ 

7.1-10“^ 

5-10“^ 

2.5 

2.54-10’ 

1.9-10“’ 

3-10'^ 

2.5 

2.54-10’ 

1.05-10“’ 

3-10“^ 

5  (L=  10  mm) 

2.54-10’ 

2.4-10“’ 

3-10“^ 

10  (L=  20  mm) 

2.54-10’ 

6.8-10“’ 

3-10"^ 

20  (L=  40  mm) 

2.54-10’ 

8.7-10“’ 

During  calculations  of  temperature  field  in  the  stationary  Stefan  problem  for 
every  fiber  diameter  the  corresponding  melting  temperature  was  calculated,  and  the 
temperatures  at  boundaries  of  liquid  and  solid  phases  were  found  by  adding  or  by 
subtracting  of  some  overheating  and  overcooling  temperature,  for  example,  100 
and  200  °C. 
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L  =25  MM,  d  =  I  |im 


U 

o 


X,  m 


Fig.  6.  Temperature  distribution  at  J  =  1  pm,  L  =  25  mm.  Pulling  velicities  1  pm/s, 
0.01  mm/s,  0.1  mm/s,  0,25  mm/s,  1  mm/s. 


Fig.  7.  Temperature  distribution  at  <7  =  1  pm,  L=  10  mm.  Pulling  velicities  1  pm  /s, 
0.01  mm/s,  0.1  mm/s,  0.25  mm/s,  1  mm/s. 
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L  =  10  mm,  d  =  0.5  |im 


U 

o 


2400  r 
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X,  m 


Fig.  8.  Temperature  distribution  at  6?  =  0.5  pm,  L  =  10  mm.  Pulling  velocities  1 
pm/s,  0.01  mm/s,  0.1  mm/s,  0.25  mm/s,  1  mm/s. 


Fig.  9.  Temperature  distribution  at  (i  =  2  pm,  L  =  10  mm.  Pulling  velosities  1  pm/s, 
0.01  mm/s,  0.1  mm/s,  0.25  mm/s,  1  mm/s. 
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2.  Determination  of  the  melting  point  of  the  composite  with  a  model. 


This  section  provides  the  basic  principles  of  the  model  of  melting  point  of 
eutectic  solid  phase.  The  basic  principle  of  the  model  for  determining  the  melting 
point  of  eutectic  is  the  concept  of  internal  energy  sources:  during  heating 
energy  of  binding  between  the  components  of  the  composite  is  released  (removed), 
as  well  as  the  internal  energy  of  thermal  stresses.  It  is  spent  for  internal  heating  of 
the  composite  (i.e.  its  self-heating)  that  causes  lowering  of  the  melting  point  of  the 
composite  (compared  to  the  melting  temperatures  of  individual  components). 

Determining  of  the  temperature  of  the  solid  phase  is  carried  out  after  the 
formation  of  the  composite  structure  (see  Fig.  10).  Therefore  depends  on  the 
composition,  structure  and  properties  of  the  composite  components  and  their 
binding  energy.  Further,  melting  temperature  will  be  taken  equal  to  solidification 
temperature,  i.  e.  Tc  =  7^.  In  reality,  it  is  an  assumption:  the  boundary  between  the 
phases  in  the  phase  diagram  (see  Fig.  10)  has  a  thickness. 


1 


T 


0 


Phc.  10.  Schematic  phase  diagram  of  two-phase  eutectic.  -  crystallization 
temperature  (solidification),  -  eutectic  concentration  of  the  phase  B ,  T^=Tj^  , 
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Fig.  11  provides  a  mechanical  interpretation  of  the  concept  of  internal 


energy  sources. 


A 


Fig.  11.  Mechanical  interpretation  of  the  concept  of  internal  energy  sources. 

The  following  spatial  energy  density  and  internal  energy  sources  are 
possible:  3Z),  ID,  D,  OZ)  (volume,  surface,  linear,  point),  see  Fig.  12. 


By  the  time  the  simplest  source  is  an  instantaneous  one.  It  is  well  known  that 
the  instantaneous  point  source  increases  the  temperature  at  the  point  to  infinity. 
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Therefore,  in  the  vieinity  of  this  souree  there  is  always  a  zone  of  melting,  the  size 
of  whieh  is  determined  by  its  power  and  ehanges  over  time. 

If  y  -  surfaee  energy  of  flat  interfaee  v  =  0  of  components  A  and  B,  then 
P(v,t)  =  y5(v)5(t)  -  instantaneous  power  of  the  source.  Depending  on  the 
conditions,  the  source  can  be  not  the  whole  boundary 
(/>(x,f)  =  y5(7-(x,f)-7;)5(f)). 


Phc.  13.  Schematic  temperature  distribution  r(v,t)  from  instantaneous  point  heat 
source  with  power  P[x,t)  =  05(v)5(t) 

In  the  general  case  source  power  depends  on  the  temperature:  P  =  P(x,t,T) . 

Approximate  determination  of  the  melting  temperature  is  made  on  the 
example  of  energy  balance  in  unit  volume  of  a  two-phase  lamellar  eutectic 
composition: 

^aPa  ~  ^o)  ~  ^aPa  (l  ~  ^0  )](^M  ~  ^0  )  ’*’  2yA^  +  W,  (24) 

where  N  =  —  -  number  of  plates  of  phase  B  with  thickness  d  per  unit  volume;  y- 
d 

surface  energy  of  the  boundary  between  the  phases  A  and  B  (it  is  assumed  that  it  is 
abruptly  released  at  temperature  T  =  T^<Tj^,  see  Fig.  14);  p^,  -  specific  density 
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of,  respectively,  phases  A  and  B\  c^,Cj^  -  specific  heat  capacity  of  phases  A  and  B\ 
q^,qjj  -  specific  latent  heat  of  fusion  of  phases  A  and  Tq  -  initial  temperature  of 
the  composite;  w  -  energy  per  unit  volume  associated  with  different  coefficients  of 
thermal  expansion  of  components: 

E,E,E,  £.  =  (l-c„)£,+c„£„ 

^  =  E^  +  \{l-CQ)aJ,E^  +  c^a^E^'\^T,  AT  =  T^-T^,  Aa  =  a^-aj, 
if  approximately  assume  that  it  is  released  at  the  temperature  T  =  .  In  (1.25) 

E^,Ej^  -  elastic  Young  moduli  of  the  composite  components.  Approximate 
estimate  (25)  for  w  is  derived  from  relations  of  linear  thermoelasticity. 
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Fig.  14.  Schematic  dependence  of  y  from  temperature. 
From  (24)  we  have  a  melting  temperature  of  eutectic: 


T  =T  =T  + 


(1b9bCo  +  (1  -  Co  )  -  ^ 

_ a 

CbPbCq  ■*"  C^Pa  (l  ~  Cq) 


(26) 


or  (when  =  p ): 


rj.  ^5Co  +  ^^(l-Co)-2yCo/^/-w  q^  q^ 

^o)I^A  ^A  ^B 


(27) 


where  and  7g  -  melting  temperatures  of,  respectively,  pure  components  A  and 

B.  From  formulas  (26),  (27)  it  follows  that  the  melting  temperature  T,n  depends  on 
the  composition,  structure  and  properties  of  composite  components  and  their 
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binding  energy.  A  comparison  of  formula  (27)  with  the  formula  of  Bystrenko  A.V. 
for  the  temperature  of  crystallization  Tc  =  Tm  implies  the  following  interesting 
relationship 


Qcq(1- Cq)  =  2ycQ/(i  +  w,  (28) 

where  Q  -  energy  of  component  mixing  in  the  melt. 

Note.  It  is  easy  to  see  that  in  determining  the  temperature  the  eutecticity 
of  the  composition  is  not  used  and  all  the  principles  of  the  model  can  be  applied  to 
any  composition,  not  only  eutectic  (i.e.  instead  of  eutectic  concentration  cq  any 
concentration  cg  [0;1]  can  be  considered).  Therefore,  with  the  proposed  model  the 
eutectic  melting  temperature  can  be  defined  as 

where  Cq  -  eutectic  concentration. 

In  the  last  formula  functional  dependencies  of  T  (□)  from  (26),  (27)  can  be 

used,  or  refined  definition  of  melting  temperature  from  the  solution  of  the 
corresponding  initial-boundary  problems  on  the  cell,  which  represent  a  composite. 

The  presented  concept  of  internal  energy  sources  is  the  basis  for  the 
formulation  of  initial-value  problem  for  the  temperature  distribution  in  the 
cell,  which  represents  eutectic  composition. 

This  problem  will  determine  the  mechanisms  of  melting  temperature 
formation  and  the  effect  on  it  of  the  composite  structure  (including  eutectic 
phase  concentration,  which  will  help  link  such  models  with  models  of  phase 
fields,  such,  for  example,  as  relation  (28)). 

The  simplest  initial-boundary  value  problem  for  monomaterial  that  explains 
the  concept  of  internal  energy  sources  (in  this  case,  source  is  in  the  material  B  and 
on  its  boundary  x  =  0,  see  Fig.  15)  and  the  nature  of  formation  of  melting 
temperature  of  eutectic: 
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(c(r)  +  ,r8(r-r„"))p 


+Y6(()8(x)  +  H<5(f), 

uX  J 


dt  dx 


V 


(29) 


XG  [0;oo), t>0 


atr<r^, 
at  T>T^  ’ 


atr<rj, 

atr>r; 


P  P5  Pl  • 

Initial  condition  for  (29):  r(a:,0)=r^,  0<a:<oo 

Boundary  conditions  for  (29):  — 

dx 

Notations: 

p  =  P5=p^  -  specific  density  of  material  (for  simplifieation  assumed  that  densities 
of  liquid  and  solid  phases  are  equal); 
c(r)  -  speeifie  heat  eapacity  of  material; 

^(r)  -  heat  eonduetivity  of  material; 
q  -  speeifie  latent  heat  of  melting; 

y5(t)5(a:)  -  density  of  heat  souree  from  binding  of  material  (e.  g.  with  another 
material)  at  the  boundary  a:  =  0,  y  =  const ; 

w5(t)  -  density  of  heat  sourees  from  internal  thermal  stresses  in  material; 

-  temperature  when  binding  at  the  boundary  a:  =  0  breaks  down  and  internal 

thermal  stresses  disappear  (see  Fig.  14); 

5(D)  -  delta  funetion; 

-  melting  temperature  of  material  B  (see  Fig.  15). 

Figure  15  gives  sehematie  graph  of  solution  T{t,x)  of  the  problem  (29)  for 

some  time  moment  t  >  0 .  It  shows  that  in  the  neighbourhood  of  the  boundary  v  =  0 
(i.  e.  in  the  neighbourhood  of  the  energy  souree)  melting  of  the  material  takes  plaee 
with  gradual  deerease  of  temperature. 


0,  T(^,t)  =  T^ 


x=0.f>0 
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Fig.  16  presents  schematic  graph  of  temperature  distribution  in  the 
neighbourhood  of  interface  between  phases  at  some  time  moment  t  >  0  after 
binding  between  phases  A  and  B  breaks  down  at  the  temperature  .  This  problem 

consists  from  two  problems  of  (29)  type. 


Note  (to  Fig.  16).  The  source  at  the  interphase  boundary  has  power 
=  0o^(O^(-^)  •  Then  in  the  neighbourhood  of  boundary  x  =  0  we  have  an 
approximation: 


Qa 


CaP. 


J 


na\t 


exp 


.2  A 


Aaltj 


Qb 


CbP, 


.4 


nalt 


exp 


,2  y 


Ao-lt  j 


«A=- 


CaPa 


2 

-  Cl  —  — — 


(30) 


^bP 


B 


Qo-Qa"^  Qb 


(31) 
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Fig.  16.  Schematic  graph  of  themperature  distribution  in  the  neighbourhood  of  flat 
boundary  x  =  0  of  initial  separation  of  phases  A  and  B  ,  (0)  =  (0)  =  0 ; 

^  _ _ V^aPa^a _ .  ^  _ _ sI^bPb^b _ 

V^aPa^a  'J^bPb^b  V^aPa^a  V^bPb^b 


Under  condition  (0,t)  =  (0,t) ,  temperature  continuity  at  the  boundary  x  =  0, 

we  have: 

Qa_Qb 
^aP  ^bPb^J^B^ 

From  this  and  from  (1.31)  we  obtain: 
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a 


Qo  =  ^aQo’  Q. 


V^bPb^Z 


V^aPa^A 


■Qo  ~  ^B0O  • 


Conclusions 

The  main  outcome  ("know-how")  of  this  section  is  the  concept  (principle)  of 
internal  energy  sources,  which  are  located  at  the  boundaries  and  in  the  bulk  of 
composite  phases  and  is  the  energy  of  binding  between  the  components  of  the 
composite  in  the  solid  state:  surface  energy  of  interphase  boundaries,  point  defects, 
internal  stresses  caused  by  the  difference  of  CTE  of  the  composite  components  and 
so  on.  The  density  of  internal  energy  sources  (it  depends  on  the  structure  of  the 
composite  in  solid  state)  may  be  such  as  to  cause  self-heating  of  composite  and 
reduce  energy  required  to  melt  the  entire  composite  when  it  is  heated,  i.  e.  reducing 
the  overall  temperature  to  melt  the  composite,  as  is  the  case  in  eutectic 
compositions.  This  reduction  can  be  significant  and  it  depends  on  the  ratio  between 
concentrations  of  the  composite  components. 

In  the  section,  we  provided  first  approximate  estimates  of  the  effect,  and 
formulated  simple  initial-boundary  value  problems  that  need  to  have  numerical 
methods  for  their  solution  be  developed  in  the  near  future  as  well  as  appropriate 
software.  For  a  more  accurate  and  thorough  investigation  of  these  effects,  a 
formulation  on  these  principles  (concept  of  internal  energy  sources)  of  initial¬ 
boundary  value  problems  for  a  representative  cell  of  the  composite  should  be 
carried  out  (some  of  such  formulations  have  been  done  by  the  authors  of  the 
report).  Special  attention  in  this  work  should  be  given  to  such  a  test  setting,  for 
which  experimental  control  will  be  possible  in  time  for  the  heating  of  the  sample 
and  temperature  distribution  in  it,  and  the  process  of  evolution  of  interphase 
boundaries  during  heating  (see  next  section). 
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Calculations 

The  estimate  of  melting  temperature  of  the  composite  by  formula  (27): 

T  =T  =T  H - - - 

^aPa  (1  ^o) 

where 


(L_^o)£<A^^  2^j.2 
2P" 

P  =  +  [(1  -  ^^0  )  ^B^A  +  C^^A^b  ]  , 


=(l-Co)^A+Co^5’ 

AT  =  T^-Tq^,  Aa  =  a^-aj, 


Table  4.  Dependence  of  calculated  melting  temperature  from  fiber  diameter 
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Fig.  17.  Calculated  dependence  of  melting  temperature  from  fiber  diameter  for 
LaB6-ZrB2  composite. 
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Fig.  18.  Enlarged  fragment  of  Fig.  17.  Dependenee  of  melting  temperature  from 
fiber  diameter  for  LaB6-ZrB2  eomposite. 

3.  Modeling  of  evolution  of  interphase  boundary  during  heating. 

In  [12]  the  results  of  an  experimental  study  were  analyzed  of  the  evolution  of 
interphase  boundaries  of  simplest  binary  euteetie  compositions  in  the  transition 
from  solid  to  liquid  state  when  heated  to  the  melting  point.  The  analysis  showed 
that  evolution  is  going  from  the  clear  boundaries  of  the  solid  phase  (i.  e.  the  surface 
between  the  phases),  to  blured  boundaries  (i.  e.  to  the  volume  boundaries  that  have 
thickness  and  are  formed  from  a  liquid  that  is  a  mixture  of  composite  components, 
see  Fig.  19  -  21). 

This  is  an  experimental  proof  of  concept  of  internal  energy  sources  according 
to  which  melting  begins  with  the  destruction  of  boundaries  between  the  phases  in 
the  solid  state.  The  controlling  mechanism  of  that  evolution  is  the  interdiffusion  of 
composite  components  (possibly  termodiffusion: 

div  p  Dgrad  c  H - grad  T  =  0 ),  which  can  form  an  eutectic  composition. 

_  V  ^  )_ 
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Next  are  shown  electron  micrographs  from  the  above  work  [12]  (Fig.  19-21) 
of  diffuse  and  eutectic  melt  zones  formed  during  heating  of  studied  composites  in 
solid  state.  The  explanations  are  provided  below  the  figures.  Detailed  explanations 
can  be  found  in  that  work. 


Fig.  19.  Electron  micrograph  of  diffuse  and  liquid  eutectic  zones  formed  by  the 
interaction  of  silver  particles  and  amorphous  silicon  film  during  annealing  {T  = 
850  °C)  [12]. 


(lOO)Si 

20nni 


Fig.  20.  HREM  micrograph  of  oriented  formation  of  eutectic  melt  in  the  system 
silver  particle  -  crystal  (100)  of  silicon  {T  =  850  °C)  [12]. 
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On  the  basis  of  analysis  of  the  results  of  direet  experimental  observation  of 
the  evolution  of  the  interphase  boundaries  of  simplest  euteetie  (and  non-euteetie) 
binary  eompositions  during  transformation  from  the  solid  state  into  a  liquid  when 
heated  [12]  and  the  eoncept  of  internal  energy  sourees  a  simplified  model  is 
eonstrueted  of  the  transition  from  surfaee  boundary  between  solid  phases  to  the 
volume  one  whieh  have  solid  and  liquid  parts. 


Fig.  21.  Eleetron  mierograph  of  diffuse  zone  and  euteetie  alloy  formed  by  the 
interaetion  of  eopper  partieles  with  amorphous  silieon  foil  under  isothermal 
eonditions  (r=810°C)  [12]. 

Main  assumption:  dijfusion  does  not  affect  temperature  distribution,  but  it 
strongly  depends  on  the  temperature. 

Therefore,  further  we  suppose  that  the  temperature  distribution  T[x,t),t>0 
in  the  vieinity  of  the  plane  interfaee  is  known  (see  Fig.  18).  Also  known  are 
dependeneies  (t),  Xj^(t),  t>0.  These  values  are  determined  by  the  solutions  of 

Stefan  type. 

Model  of  phase  interdiffusion. 

It  is  assumed  that  interdiffusion  is  determined  (mainly)  by  the  dependenee  of 
the  diffusion  eoeffieient  on  temperature 
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D,{T)  =  D,.&xp 


r(t-t^) 


,  i  —  A,  B  y 


(32) 


where  subseript  i  indieates  phase,  whieh  diffusion  is  eonsidered. 

Instead  of  dependence  (32)  any  other  dependence  can  be  used.  Examples  of 
direct  experimental  observation  of  processes  of  diffusion  and  evolution  of 
interphase  boundaries  of  simplest  eutectic  (and  non-eutectic)  binary  compositions 
during  the  transition  from  solid  state  to  liquid  when  heated  can  be  found  in  [12]. 

Approximately,  we  assume  that  in  the  area  of  the  melt  <  v< 

phases  are  distributed  homogeneously  due  to  high  values  of  diffusion  coefficients 
at  high  temperature  (see  Fig.  22,  where  a  scheme  of  distributions  of  phase 
concentrations  in  the  vicinity  of  the  interface  is  given).  This  assumption  can  be 
weakened  by  the  assumption  of  their  linear  distribution  (with  a  constant 
concentration  gradient  other  than  0).  However,  it  should  be  noted  that  this 
approach  to  the  study  of  diffusion  in  the  melt  is  very  approximate.  Diffusion  in  this 
area  is  a  complex  process,  which  is  linked  to  other  processes,  and  requires  more  in- 
depth  research  and  modeling. 


Fig.  22.  Schematic  distribution  of  concentration  of  phases  A  and  B  for  some  time 
moment  t>0.  l{t)  -  length  of  melt  zone. 
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Further  diffusion  problem  (for  shortness)  will  be  formulated  for  phase  A.  For 
phase  B  the  problem  is  formulated  similarly.  Thus,  in  the  region 
(t)  < X <  (t))  for  eoneentration  C^[x,t)  we  have 


C,(x,t)  =  C,(x„t)  =  C,(Oj)  =  C,(x„t)  =  C,(t),  C,(0)  =  1.  (33) 

In  the  region  (t) < v < (t))  eoneentration  C^[x,t)  is  found  from 
equation 


dx 


dx 


(34) 


aeeording  to  whieh  eoneentration  field  C^(x,t)  is  eonsidered  in  quasistatie 
approximation.  Right  end  of  region  <  v< is  defined  by 

temperature  ,  at  whieh  diffusion  of  A  eomponent  into  B  eomponent  stops,  i.  e. 
from  equation 

T{x^aA)  =  Tj,^-  (35) 


From  (32)  for  density  of  diffusion  flux  F(t)  of  phase  A  in  the  analysed  region  we 
have 


F{t)  =  -D^(T{x,t))-^.  (36) 

Integrating  the  last  equality  in  v,  for  eoneentration  C^(x,t),  (t)  <  v  <  (t) ,  we 

have 


X 

C^(x,t)  =  -F(t)  J 


dx 


.DA(T(x,t)) 


+  G{t),  t>0, 


(37) 
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where  function  G{t)  is  found  from  the  condition 


0  =  -F(t)  f  — x'  +  P(0 
Therefore,  equality  (37)  can  be  rewritten  as 


(38) 


^DAV 

C,{x,t)  =  F{t)  J 


dx 


Xg(t)<x<XgX‘)’  '-0-  (39) 


D,(T{x.t))' 

Finally,  the  distribution  C^(x,t)  will  be  known  if  we  determine  density  of  flux 

F(t)  of  phase  A.  This  density  is  determined  from  the  condition  of  conservation  of 
quantity  of  phase  A: 

^DA(f) 

l(t)-Xj^(t)  =  C^[xj^(t),t)l(t)+  j  C^(x,t)dx, 


^b(0 


which,  taking  into  account  (39),  can  be  written  as 

^  ('f'l 

dx 


i  D^(T(x,t))  J  D,{T(x,t)} 

whence  we  obtain  for  the  flux  density  F(t) 


F(t)  = 


l(t)  I  - - yr+  \  dx  {  - - y- 

'  '  i  D,(T(y,t))  J  J  D,(T(y,t)) 


^daV)  ^DA\n 


dy 


(40) 


Thus,  from  the  formulas  (39),  (40)  for  the  distribution  of  concentration  of  phase  A 
we  obtain  the  following  explicit  formula 


C^{x,t)=- 


U{t)-xJtj\ 
LU  .WJ  J 


)) 


%aW  ,  %a(^)  %a(0  , 

,(,)  f  — r  dx  I 

i  D,(T{y,t))  {  D,(T(y,t)) 


Xs{t)<x<x^^{t),  t>0  (41) 
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The  peculiarity  of  the  solution  (41)  is  its  uncertainty  at  t  =  0:  C^(0,0) 


— ,  which 
0 


strongly  depends  on  the  functions  x^{t)  defined  by  the  solution  of  the 

problem  of  Stefan  type,  as  well  as  the  solution  Xjj^{t)  of  equation  (35),  which  is 

determined  by  the  type  of  dependence  (32).  This  uncertainty  can  be  resolved  and 
examined  using  the  solution  of  the  discussed  temperature  problem  or  through 
computational  experiments. 


Fig.  23.  Simplified  scheme  of  evolution  of  flat  boundary  v  =  0  between  solid 

phases  A  and  B  (T  >T^,T^<T^ ). 

When  x>Xj^^(t)  A  phase  is  absent,  and  C^(x,t)  =  0  so  region  x>Xj^^(t) 
consists  only  of  phase  B.  In  the  case  when  A  phase  does  not  diffuse  into  solid  phase 
B,  formula  (1.40)  determines  concentration  C^(x,t)  in  the  melt: 


C^{x,t) 


<v<Vg(t),  l^(t)  =  -x^(t) .  This  case  is  shown  in  Fig.  23. 


Fig.  23  gives  simplified  scheme  of  evolution  of  planar  interphase  boundary  in 
the  case  of  coincidence  of  diffusion  and  melt  zones  (or  in  the  case  when  mismatch 
can  be  neglected)  and  notations  used:  l(t)  -  thickness  at  time  moment  t  of  the  liquid 
phase,  which  is  located  between  solid  phases  A  and  B  and  consists  of  uniformly 
dissolved  components  A  and  B.  The  assumption  of  uniformity  is  justified  in  the 
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case  of  a  higher  rate  of  diffusion  in  the  liquid  than  in  the  solid  phase.  (t) ,  4  (t)  - 
thicknesses  of  phases  A  and  B  dissolved  in  the  melt  at  time  moment  t,  respectively. 

In  this  case,  we  can  assume  that  Q  =  y  and  =  (1  - Q)  =  y  are  concentrations 


in  the  liquid  of  phases  A  and  B,  respectively. 

Thus,  main  formulas  of  the  model  are  as  follows. 

For  the  distribution  of  concentrations  C^[x,t),Cg[x,t)  of  phases  A  and  B 
under  heating  we  have: 


C^ix,t)-- 


C^(xj^{t),t),  0<x<Xi^{t),  t>0; 

^da(0  r 

[/(/)-Ag(/)]  f  — ^ ^ 
i  D,{T{y,t)) 


ml 


^da(0  r  ^daW  ^daW  r 

_ ^ _ +  r  dx  f  _ ^ _ 

DAHyA)  i,i  i  DAT(y,t)) 


;  Xj^{t)<x<Xj^^{t),  t>0;  (42) 


0,  x>Xj^^{t),  t>0, 


where  functions  Xg(t),l(t),T{^x,t)  -  are  determinedfrom  corresponding 

temperature  problem  of  Stefan  type,  and  implicit  function  is  a  solution  of 

equation 

T{x^aA)  =  Tj,^-  (43) 


It  is  assumed  that  in  (32)  dependence  of  diffusion  coefficient  D^{T)  from 
temperature  is  known. 

As  was  noted  earlier,  the  feature  of  the  solution  (42)  (43)  is  its  uncertainty  at  t 


=  0:  Cg(0,0) 


— ,  which  strongly  depends  on  the  functions  Xsit),  l(t),  T(x,t),  XoAit). 


This  uncertainty  can  be  resolved  and  investigated  with  the  help  of  the  previously 
mentioned  temperature  problem  of  Stefan  type  or  computational  experiments. 
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Conclusions 


This  section  presents  a  new  model  of  the  evolution  of  the  interphase  boundary 
during  heating,  according  to  which  surface  boundary  between  the  phases  is 
transformed  into  a  boundary  that  has  thickness  (is  volumetric)  and  consists  of  two 
solid  parts  between  which  is  liquid.  The  model  was  developed  on  the  assumption 
that  diffusion  does  not  affect  temperature  distribution  which  is  determined  through 
solution  of  separate  problem,  but  temperature  affects  diffusion  through  diffusion 
coefficient  dependence  on  temperature. 


4.  Phase-field  study  of  structure  formation  in  binary  eutectic 

The  goal  of  this  work  is  the  theoretical  study  and  computer  simulation  of  the 
processes  of  structure  formation  and  the  properties  of  interfaces  in  binary  eutectic 
systems.  As  the  result  of  the  thermodynamical  instability  of  the  uniform 
composition,  these  systems  tend  to  form  spatial  structures  in  the  solid  state,  the 
eutectic  colonies  [13-15].  The  significance  of  the  study  of  the  properties  of  eutectic 
systems  for  material  science  is  due  to  wide  industrial  applications  of  these  latter. 
This  concerns,  in  particular,  the  processes  of  sintering  in  powder  metallurgy,  the 
utilization  of  the  directional  crystallization  in  manufacturing  fiber  composites  for 
electronics  and  photonics,  etc.  The  investigation  of  the  properties  of  interfaces  and 
the  closely  related  phenomenon  of  contact  melting  are  therewith  the  central  points 
for  understanding  the  processes  of  structure  formation  in  eutectic  systems. 

Within  the  framework  of  the  project,  the  following  issues  were  considered: 

(1)  formation  of  eutectic  structures  from  the  overcooled  melt  (glass)  obtained  by 
the  rapid  quenching  of  the  liquid  LaBg  -ZrB2  eutectic; 

(2)  formation  of  eutectic  colonies  in  the  process  of  directional  crystallization  of 
LaBg  -  ZrB2  ceramics; 
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(3)  study  of  the  fundamental  interfaee  properties  and  the  phenomenon  of  eontaet 
melting  in  binary  euteetie  taking  as  an  example  the  Cu-Si  system. 

As  the  investigation  tool,  the  phase  field  theory  (PFT)  [16-18]  with  the 
dissipative  mierostruetural  alloy  dynamies  deseribed  by  the  parabolie  equations  for 
phase  variables  determining  the  loeal  phase  (liquid  or  solid)  system  state  has  been 
employed.  The  basie  idea  of  PFT  is  to  eonsider  the  inter-phase  boundaries  as  the 
eontinuous  spatial  regions  of  finite  width,  the  approaeh,  whieh  essentially 
simplifies  the  numerieal  ealeulations. 


PFT  equations  and  the  model  for  binary  eutectic 

Within  the  PFT,  the  dissipative  kineties  of  a  system  is  deseribed  in  terms  of 
phase-field  variables.  The  relevant  kinetie  equations  for  an  isothermal  binary 
system  ean  be  obtained  from  the  free  energy  functional. 


f{cp,c,T)^^{Vcf^^{Vcpf 


dV , 


(44) 


where  f{(p,c,T)  is  the  free  energy  density;  T,  c,  and  (p  denote  the  temperature, 
concentration  and  the  phase  field  variable,  respectively.  The  phase-field  variable 
specifies  the  local  state,  ^  =  0  for  solid,  and  (p=\  for  liquid  phase.  By  applying  the 
basic  principles  of  thermodynamics,  one  obtains  the  PFT-equations  (with 
allowance  for  the  convective  medium  motion  with  the  velocity  Vg  along  the  OX 
axis),  in  the  form  [16-17] 


dcp 


-M. 


+  v, 


dcp 

dx 


(45) 


^=V- 

dt 


MpC(l-c)V 


ac  " 


+  Vn 


dc 

dx 


(46) 


yoc  j_ 

Here  the  quantities  and  are  the  kinetic  coefficients  specifying  the 
relaxation  rate  of  the  system  to  the  equilibrium;  ande^  are  the  parameters 
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responsible  for  the  surfaee  energy  of  solid-liquid  and  solid-solid  interfaees.  The 
speeifie  features  of  a  system  are  determined  by  the  partieular  form  of  free  energy 
density /(^9,c,r).  In  this  work  we  use  for  the  latter  the  standard  expression  for  a 
non-ideal  binary  system  with  the  miseibility  gap  [16] 

c,  T)  =  (l-  c)Ug>,  T)  +  cf,  (g>,T) 

RT  ( 471 

+ - [clnc  +  f  1  -  cjlnf  1  -  cj]  +  cf  1  -  c)[Qs  [1  -  p{  f)]  + 

The  typieal  minimum  in  the  melting  eurve  on  the  phase  diagram  is  therefore  due  to 
the  differenee  in  mixing  energies  for  liquid  and  solid  phases.  Here  R  is  the  gas 
eonstant,  is  the  molar  volume,  and  are  the  mixing  energies  for  the  solid 
and  liquid  state,  respeetively.  The  eomponent  free  energies  and  fei^pj) 

have  the  form 

/„(«.,r)  =  (48) 

where  T’m  *  are  the  melting  temperatures,  are  the  energy  barriers  assoeiated 
with  the  surfaee  energy  of  liquid-solid  interfaees,  are  the  eomponent  latent 
heats  for  the  eomponents  A  and  B;  the  funetions  g{(p)=  (p^{l-(pf  and 
p{(p)=  ^9^(l0-15^9  +  6(p^)  are  the  model  barrier  and  interpolating  funetions  eonstrueted 
in  sueh  a  way  as  to  deseribe  the  liquid-solid  interfaees  of  a  finite  width.  The  kinetie 
eoeffieients  and  are  determined  for  a  binary  system  as 


M^  =  (\-c}M^  +  cM, 


M 


c  ~ 


Ds  +  P((P)(Dl-Ds) 
RT/v^ 


(49) 


with  Dg  and  being  the  solid  and  liquid  diffusivities,  respeetively. 

As  was  shown  in  reeent  experiments  with  formation  of  diffusion  zone  and 
eontaet  melting  in  binary  euteetie,  the  equilibrium  interfaee  ean  exhibit 
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complicated  (three-phase)  structure  and  the  width  of  the  order  of  several  micron 
[19,  20].  This  means  that  the  interface  width  (as  well  as  the  surface  energy)  is  a 
real  physical  quantity  (but  not  the  model  parameter,  as,  for  instance,  in  Ref.  [18]), 
which  has  to  be  determined  from  the  experiment. 

For  numerical  solution  of  the  equations  (45-46)  there  has  been  employed  in 
all  cases  the  finite  element  method  implemented  in  the  free  program  package 
FEniCS  [21]. 

Simulation  of  the  crystallization  of  eutectic  structures  from  the  overcooled  melt 

To  investigate  the  crystallization  of  the  eutectic  from  the  melt  obtained  as 
the  result  of  rapid  quenching,  the  PFT  equations  (45-46)  have  been  solved  within 
the  isothermal  approximation.  The  initial  state  of  a  binary  system  has  been 
specified  as  the  overcooled  liquid  state  (i.e.  as  the  eutectic  glass);  the  composition 
was  set  uniform,  with  constant  concentration  and  the  Gaussion  noise  of  0.025- 
0.035  amplitude,  with  the  distribution  of  the  phase  variable  corresponding  to  the 
liquid  {(p  =  l)  or  solid  {(p  =  0)  state.  The  temperature  was  set  below  the  eutectic 
temperature  in  the  area  of  the  solid  state  with  the  eutectic  composition;  the 
homogeneous  Neumann's  boundary  conditions  were  used. 

The  model  PFT  parameters  were  chosen  to  be  as  close  as  possible  to  the 
available  experimental  data  for  the  real  binary  LaBg-ZrB2  ceramics  [15],  or  to 

reproduce  the  known  component  properties  (Table  5).  For  instance,  the  latent  heats 
and  mixing  energy,  which  determine  the  equilibrium  properties  of  the  system,  were 
adjusted  to  reproduce  the  eutectic  point  on  the  phase  diagram;  the  kinetic 
coefficients  were  estimated  from  the  available  experimental  data  on  the  velocity  of 
the  directional  solidification;  the  values  for  the  gradient  coefficients  were  chosen  to 
reproduce  the  typical  geometrical  parameters  of  eutectic  colonies. 
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Table  5.  Parameters  for  the  PFT  model  used  in  simulations  of  euteetie 


erystallization  from  the  overeooled  LaBg  -  ZrB2  melt. 


Parameter 

Value 

Units 

Density,  p 

LaBg 

4730 

kg/m^ 

ZrB2 

6090 

Simulation  eell  size,  X 

0.5 -40.0  •10“'’ 

m 

Melting  temperature, 

LaBe 

3000 

°K 

ZrB2 

3310 

Kinetie  eoeffieient  for  the  inter¬ 
phase  boundary,  =  pg 

1,0- 10'’ 

m  /(°K  see) 

Liquid  and  solid  diffusivities. 

1,0-10-'",  MO''^ 

m^/see 

Surfaee  energy  for  the  liquid- 
solid  interfaee,  =  Og 

5,0 

J/m" 

Molar  volume 

LaBg 

43,1-10'’ 

m^/mole 

ZrB2 

16,35-10'’ 

Euteetie 

35,3-10''’ 

Euteetie  temperature  for 

LaBg  -  ZrB2 

2740 

OR 

Latent  heats,  L 

LaBe 

3200-10'’ 

J/m^ 

ZrB2 

7200-10'’ 

Mixing  energy  for  solid  state, 

2800-10'’ 

Euteetie  eoneentration 

0,32 

- 

In  numerieal  solution  of  PFT  equations  the  spatial  domain  was  ehosen 
within  the  range  of  X=0.5-40.0  mieron;  the  liquid-solid  interfaee  width  (the  model 
parameter)  was  set  5/X  =0.05.  The  absolute  error  of  eomputations  for  the 
eomposition  and  phase  profiles  was  kept  less  than  10'^ . 

All  the  parameters  entering  PFT  equations,  whieh  are  neeessary  for  solving 
the  problem,  ean  be  expressed  via  the  quantities  given  in  Table  5  by  using  the 
known  relations  [16]. 
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Fig.  24.  The  results  of  ID-simulations  of  eutectic  structure  formation  from  glassy 
state  in  LaBg  -ZrBj  for  the  temperature  T  =  2000° for  the  time  points  t=0;  4;  280 
sec  (from  the  left  to  the  right). 


Fig.25.  The  results  of  2D- simulations  (composition  distribution)  of  the  formation 
of  eutectic  structure  in  LaBg  -  ZrBj  from  glass  at  the  temperature r  =  2000° for  the 
time  points  t=0;  60;  180  sec  (from  the  left  to  the  right). 

The  results  of  ID-  and  2D-simulations  of  the  behavior  of  the  system  with  the 
eutectic  composition  after  the  rapid  quenching  to  the  temperature  r  =  2000°^  are 
given  in  Figs.  24  and  25.  It  should  be  mentioned  that  at  the  initial  stage  the  system 
crystallizes  into  the  weakly  structured  solid  phase,  therefore,  the  process  of  the 
formation  of  colonies  occurs  in  the  solid  phase,  that  is,  actually,  from  glass. 


Distribution  A:  Approved  for  pubiic  reiease;  distribution  is  uniimited. 


The  above  simulations  reproduce  the  process  of  the  formation  of  the  spatial 
structures  as  a  whole.  It  is  worth  mentioning  that  the  series  of  very  long  runs 
(f>10^  sec)  have  demonstrated  the  capability  of  the  model  to  reproduce  the 
complete  decomposition  of  the  system  into  components  while  approaching  the 
final  equilibrium  state. 

Simulation  of  the  directional  crystallization  of  fiber  structures  in  ceramics 

The  directional  crystallization  was  studied  in  the  "frozen  temperature 
approximation",  i.e.  the  temperature  distribution  was  supposed  to  be  steady  and 
linear  along  the  direction  of  sample  motion.  Such  an  approximation  is  justified, 
since  the  temperature  conductivity  coefficient  of  boride  ceramics  considerably 
exceeds  the  mass  transport  coefficients.  The  initial  phase  distribution  was  set  in  the 
form  of  the  smoothed  Heaviside  function  describing  the  liquid-solid  interface  with 
the  front  at  the  eutectic  temperature.  The  initial  composition  distribution  was  set 
equal  to  the  constant  eutectic  concentration  with  the  Gaussian  0.025-0.035 
amplitude  noise  (Fig.26-27) 

The  simulations  have  been  performed  for  parameters  given  in  Table  5.  The 
spatial  domain  of  simulations  was  within  the  range  X=20.0-40.0  mkm;  the 
temperature  gradient  was  set  equal  to  4.5  ■  10^  K/m;  the  velocity  of  the  directional 
crystallization  was  in  the  range  =0.25 -2.0  10^^  m/sec. 

concentration 

-^0.325 

I “0,323 
0.320 
0.318 
0.315 


Fig.26.  Initial  distributions  for  the  concentration. 
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Fig.28.  The  results  of  2D-simulations  for  the  phase  (top)  and  concentration 
(bottom)  for  the  initial  stage  of  the  process  of  the  directional  crystallization  for  t=\  \ 
2;  12  sec  (from  the  left  to  the  right).  The  velocity  of  the  melt  motion  (from  the 
right  boundary)  is  v  =  0.5- 10“^  cm /sec . 
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Fig.29.  The  results  of  the  2D-simulations  of  the  directional  crystallization  for  the 
phase  (top)  and  concentration  (bottom)  for  the  time  point  t=10  sec  and  the  melt 
velocity  (from  the  right  boundary)  V  =  0.5  10"^  0.75- 10“^;  1.25- 10“^m/ sec  (from 
the  left  to  the  right). 

The  results  of  2D-simulation  of  the  directional  crystallization  in  boride- 
boride  ceramics  are  given  in  Fig.  28.  At  the  initial  stage,  there  emerges  an 
instability  on  the  crystallization  front,  which  gives  rise  to  the  composition 
segregation  and  the  formation  of  a  cell  structure  along  the  front.  Then,  at  the 
transitional  stage,  a  part  of  cells  gets  destroyed  and  the  eutectic  structure  growth 
larger.  Finally,  due  to  the  transition  processes,  there  forms  a  steady-state  fiber 
structure  representing  the  eutectic  colony  being  the  result  of  the  directional 
crystallization.  The  results  are  in  agreement  with  the  known  experiments  [15]. 
Notice  that  the  simulations  clearly  demonstrate  the  stabilizing  effect  of  the 
directional  melt  motion  on  the  structure  formed.  In  particular,  the  fact  of  the 
existence  of  the  minimal  velocity  needed  for  the  directional  crystallization  is 
reproduced. 
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Fig.30.  Dependence  of  the  space  parameter  of  the  fiber  structure  on  the  velocity  of 
the  directional  crystallization  of  LaBg  -  ZrB2  .The  results  of  PFT-simulations  are 
marked  as 

The  simulation  results  for  different  velocities  of  directional  crystallization 
are  given  in  Fig.29.  As  follows  from  the  data  obtained,  the  typical  space  parameter 
of  the  fiber  structures  formed  depends  on  the  velocity.  Growing  velocity  of  the 
directional  crystallization  results  in  smaller  colonies.  The  relevant  dependence  is 
given  in  Fig.  30.  Also,  the  simulations  reproduce  the  existence  of  the  maximal 
threshold  velocity  of  the  directional  crystallization,  above  which  the  structures  do 
not  form.  In  these  calculations  this  threshold  was  estimated  as  v  =  1.75  10^^m/sec. 
Apparently,  this  is  due  to  the  fact  that  rate  of  the  diffusive  processes  in  this  case 
turns  out  to  be  insufficient  to  provide  the  component  segregation.  Let  us  note  that 
these  results  correlate  qualitatively  with  the  known  experiments  [15]. 
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Study  of  the  contact  melting  and  the  structure  of  binary  eutectic  near  the  eutectic 
point 

The  study  of  the  interfaee  phenomena  in  binary  euteetie  was  aimed  basieally 
at  the  phenomenon  of  eontaet  melting  (CM)  and  the  properties  of  diffusion  zone. 
The  goal  was  to  obtain  the  interpretation  for  the  results  of  reeent  experiments  with 
binary  systems  eonsisting  of  metal  Au,  Al,  Ag,  and  Cu  partieles  (of  ==5  10  '’m  size) 
plaeed  on  amorphous  or  erystalline  silieon  films  [19,  20].  The  general  eonelusions 
obtained  in  these  works  ean  be  formulated  as  follows:  1)  for  the  temperatures 
below  the  euteetie  point,  under  isothermal  eonditions,  the  solid-solid  interfaee 
forms  a  diffusion  zone  (DZ)  with  typieal  width  of  ~  10^  m  ,  whieh  may  reduee  in 
size  for  lower  temperatures;  2)  for  the  temperatures  within  the  range  10-20°C 
above  the  euteetie  point,  under  isothermal  eonditions,  the  formation  of  the  liquid 
euteetie  layer  at  the  interfaee  (i.e.,  the  phenomenon  of  CM)  is  observed,  whieh  is, 
however,  always  preeeded  by  the  formation  of  the  DZ. 

In  simulations  there  has  been  used  the  standard  model  formulated  in  See.  2. 
It  should  be  noted  that  the  basie  idea  underlying  all  the  following  eonsideration  is 
that  the  DZ  ean  be  viewed  as  the  diffuse  solid-solid  interfaee  assoeiated  with  the 
state  of  the  equilibrium  with  eomplete  deeomposition.  Furthermore,  we  negleet  the 
effeets  of  a  large  number  of  possible  eompounds  in  DZ  with  different  struetural 
and  thermodynamieal  properties,  and  use  the  eoneentration  field  as  the  only 
variable  deseribing  the  eomposition. 

We  give  here  the  results  of  PFT  simulations  for  Cu-Si  system.  The 
following  parameters  were  used  (A  =  Cu;B  =  Si):  the  domain  of  solution 
=  2.0  10'^  m;  time  interval  =  4000  sec ;  melting  temperatures  Tf  =  1357  °K 
and  T’m  =1685  °K;  latent  heats  =  1.837  •lO*' //m^  and  =  4.488  10'' //m^ 

molar  volume  =  8.57  ■10'’ mV mo/e ;  surfaee  energies  for  liquid/solid  interfaees 
a  a  ^  =  03  J  Im^ ;  diffusivities  =  10'^  n3 !  sec  and  =  10  ''  I  sec ;  kinetie 
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parameters  for  interfaees  //^  =  //^  =  5  •  10'^  ml{°Ksec)  ;  mixing  energy  for  solid  phase 
i3  =  3.635-10‘’  Jlm^ . 

Most  of  the  above  parameters  are  the  experimental  data  available  for  Cu-Si 
system  or  the  typieal  numbers  for  the  relevant  quantities  (like  solid/liquid 
diffusivities)  used  in  similar  simulations  [22-24].  The  experimental  value  for  the 
important  quantity  Q  is  not  available  in  the  literature;  we  fitted  it  to  make  the 
euteetie  temperature  mateh  its  experimental  value  =1075°Z.  In  order  to  make 
eomputer  simulations  less  time  eonsuming,  the  kinetie  eoeffieients  were  set 
much  less  than  the  typical  values  used  in  similar  PFT  simulations  [24].  However,  it 
is  to  emphasize  that  the  simulations  performed  were  aimed  at  the  investigation  of 
the  steady  (equilibrium  or  metastable)  states,  and,  as  is  seen  from  Eq.  (46),  the 
kinetic  coefficients  affect  the  rate  of  the  dissipative  relaxation  of  the  system  rather 
than  the  final  steady  states.  The  important  for  the  model  quantity,  the  width  of  the 
DZ,  is  determined  by  the  concentration  gradient  coefficient  .  The  latter  was 
fitted  as  well  to  provide  the  typical  value  d  =  410  m  observed  experiments.  The 
associated  interfacial  energy  of  DZ  evaluated  in  simulations  was  on  the  order  of 
2’-500//m"  . 

The  above  parameters  ^  and  are  needed  to  calculate  the  phase  gradient 
coefficient  ,  the  barrier  energy  W^  g  ,  and  the  kinetic  parameters  M  ^  g  entering 

the  PFT  equations;  for  this  purpose,  the  standard  relations  given  in  Ref.  [16]  were 
employed.  The  model  parameter  specifying  the  liquid-solid  interface  width  was 
set  ^/Z_  =  0.01. 

The  absolute  error  for  composition  and  phase  profiles  obtained  in 
computations  was  kept  less  than  10^^ . 
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Fig.  31.  Experimental  phase  diagram  of  Cu-Si  in  the  vicinity  of  the  eutectic  point 
[22].  The  eutectic  concentration  is  0.3. 

A  part  of  the  experimental  phase  diagram  for  Cu-Si  near  the  eutectic  point 
[22]  is  displayed  in  Fig.  31  To  reproduce  the  results  of  experiments  with  the 
formation  of  the  DZ  [19],  we  performed  PET  simulations  in  ID  (planar)  geometry 
for  the  points  D1  (T^j  =  923°K  (650°C))  and  D2  =  1083°K  (810°c))  for  the  eutectic 
concentration  c=0.3.  For  this  purpose  Eqs.  (45-46)  were  solved  with  initial 
conditions  describing  a  contact  of  two  solid  components  (Fig.  32,  line  1). 


Fig.  32.  Initial  concentration  profiles  (solid  lines  1  and  2)  and  phase  distribution 
(dashed  stochastic  line)  used  in  PFT  simulations. 


Distribution  A:  Approved  for  pubiic  reiease;  distribution  is  uniimited. 


1 


c  0.8 
g 

4_f 

(0 


c 

<v 

u 

c 

o 

u 


0.6 


01 

to 

(D 

Q. 


0.4 


0.2 


0 


0 


_i _ I _ I _ I _ I 

0.2  0.4  0.6  0.8  1 

distance,  x/Xm 


Fig.  33.  Final  steady-state  profiles  for  phase  (dashed  line)  and  eomposition  (solid 
lines)  obtained  in  PFT  simulations  for  the  temperatures  =  923° K  (650 °C)  (line  1) 
and  =  413°K  (200°C)  (line  2)  for  the  time  point  t  =  (4000  see). 

Notiee,  that  the  PFT  equations  allow  overheated  metastable  solid  state 
(^9  =  0)  to  exist  infinitely  long.  For  this  reason,  to  trigger  the  kinetie  proeesses,  the 
initial  phase  distribution  was  modified  to  inelude  0.05-0.15  -  amplitude  random 
noise.  In  prineiple,  the  effeets  of  random  fluetuations  ean  be  taken  into  aeeount  by 
adding  noise  souree  to  PFT  equations  or  using  the  renormalization  approaeh  [25- 
27].  However,  we  did  not  eonsider  the  fluetuations  in  simulations,  since  these  latter 
affect  the  kinetics  rather  than  the  final  steady  states. 

The  results  of  simulations  for  the  temperatures  specified  in  experiments  with 
Cu-Si  system  [19]  are  displayed  in  Figs.  33,  34.  For  the  temperature 
=  923°K  (650°C)  below  the  eutectic  point  (Fig.  33),  the  simulations  show  that,  as 
the  result  of  the  relaxation,  there  forms  a  final  steady  state  corresponding  to  'solid 
A  -I-  diffuse  interface  -i-  solid  B'  composition  profile.  It  is  worth  mentioning  that 
using  any  of  the  initial  composition  distributions  given  in  Fig.  32  yields  identical 
final  steady-state  composition  profiles,  which  is  the  result  of  the  equilibrium  nature 


Distribution  A:  Approved  for  pubiic  reiease;  distribution  is  uniimited. 


of  this  state.  The  area  in  grey  (diffuse  interface)  in  Fig.  32  is  interpreted  as  a  DZ. 
Another  important  observation  obtained  in  PFT  simulations  is  the  reduce  in  the 
width  of  DZ  with  lowering  the  temperature. 

The  simulations  performed  for  the  temperature  =  10S3°K  (810°C)  above 
the  eutectic  point  (Fig.  34)  indicate  the  formation  of  a  liquid  area  within  the  steady 
DZ,  quite  consistent  with  experimental  observations  [19].  Notice  that  the  steady 
three-phase  (solid  A  -i-  liquid  -i-  solid  B)  profiles  can  be  viewed  as  the  state  of 
complete  decomposition  (like  those  shown  in  Fig.  33),  but  with  a  liquid  interface. 
It  is  to  mention  that  similar  simulations  performed  for  the  parameters  specific  for 
Ag-Si  system  yielded  quite  analogous  results. 

The  fundamental  question,  which  arises  in  this  context,  is  that  of  the  nature 
of  the  steady  states  observed,  i.e.,  if  they  are  metastable  or  equilibrium  ones.  To 
address  this  question,  we  performed  series  of  runs  with  various  initial  conditions  in 
the  vicinity  of  the  expected  eutectic  point  with  monitoring  the  free  energy  of  the 
system  defined  by  Eq.  (44).  Depending  on  the  initial  phase  and  concentration 
distributions,  for  the  same  temperature  and  average  concentration  (i.e.,  for  the 
same  point  on  the  phase  diagram),  the  system  may  be  "trapped"  in  a  metastable 
state,  or  proceed  to  the  thermal  equilibrium. 

An  example  of  the  behavior  of  free  energy  during  relaxation  for  the 
temperature  T^2  =  1083°.^:  (810°C)  and  eutectic  concentration c  =  0.3  (point  D2)  for  the 
three-phase  (DZ  with  liquid  interface)  and  two-phase  (DZ)  states  is  given  in  Fig. 
35  (with  free  energy  density  measured  in  units  of  A  =  RT^  /v„  =  1.317  •  10‘^  J I  ). 

The  free  energy  density  for  liquid  state  evaluated  for  the  same  point  on  the 
phase  diagram  is  much  higher,  f^l fo=  0.0705 .  Therefore,  one  must  conclude,  that 

1)  the  liquid  state  at  the  point  D2  is  metastable;  2)  thermal  equilibrium  in  the 
system  is  associated  with  the  steady  three-phase  'solid  A  -i-  liquid  -i-  solid  B'  state, 
because  it  has  the  lowest  free  energy. 
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Fig.  34.  Final  steady-state  profiles  for  phase  (dashed  line)  and  eomposition  (solid 
line)  obtained  in  PFT  simulations  for  the  temperature  r^2  =  1083°Z  (810°c)  (point 
D2). 


Fig.  35.  Average  free  energy  density  of  the  system  as  a  funetion  of  time  for  the 
temperature  7^,2  =  1083°Z  (810°C)  and  euteetie  concentration  c  =  0.3  for  the  steady 
states:  solid  A  -i-  solid  B  (complete  decomposition;  solid  line);  solid  A  -i-  liquid  -i- 
solid  B  (complete  decomposition  with  a  liquid  interface;  dashed  line). 
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Fig.  36.  Average  free  energy  density  of  the  system  vs.  temperature  for 
different  steady  states  for  the  euteetie  eoncentration  c=0.3.  The  area  in  grey  points 
to  the  temperature  interval,  where  the  three-phase  state  has  the  lowest  free  energy. 


Fig.  37.  Average  free  energy  density  of  the  system  vs.  temperature  for  different 
steady  states  for  the  eoneentration  e=0.4.  The  areas  in  light  and  dark  grey  point  to 
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the  temperature  intervals  of  thermodynamical  stability  of  three-phase  'solid  A  -i- 
liquid  -I-  solid  B'  and  two-phase  'liquid  -i-  solid  B'  states,  respectively. 


Fig.  38.  Phase  diagram  near  the  eutectic  point  obtained  from  PFT  simulations 


It  is  to  emphasize,  that  this  theoretical  (based  on  PFT)  conclusion,  while 
being  consistent  with  above  cited  experiments,  contradicts  the  experimental  phase 
diagram  (Fig.  31),  since  the  point  D2  lies  here  within  the  area  of  liquid  phase.  To 
examine  this  issue  in  more  detail,  we  performed  further  simulations  aimed  at  the 
evaluation  of  free  energy  of  possible  steady  states  in  the  vicinity  of  eutectic  point 
and  deriving  the  associated  phase  diagram.  The  results  are  given  in  Figs.  36-38. 

The  most  unexpected  finding  is  that,  according  to  the  PFT  model  used,  the 
solidus  line  T  =  turns  out  to  have  a  'fine  structure',  i.e.,  it  splits  and  expands  to 
the  relatively  narrow  (about  =  50  °C  width)  area  above  the  eutectic  temperature, 
where  the  above  mentioned  states  describing  the  complete  decomposition  with  the 
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liquid  interface  ('solid  A  +  liquid  +  solid  B')  are  thermodynamically  stable.  The 
point  D2  lies  therewith  within  this  area,  which  explains  the  experimental 
observations  of  the  DZ  with  the  layer  of  liquid  eutectic  [19]. 

In  this  context  a  few  words  should  be  said  about  the  standard  technique  of 
the  common  tangent  construction  for  obtaining  phase  diagrams.  The  latter  ignores 
the  gradient  terms  in  the  free  energy  functional  (44)  and,  therefore,  implies  the 
equilibrium  in  the  system  to  be  associated  with  the  uniform  concentration  and  the 
phase  field.  As  the  result,  the  contribution  to  the  free  energy  of  the  interfaces  is 
neglected  and  the  description  of  the  complicated  interfaces  like  given  in  Fig.  34 
with  high  interfacial  energy  is  impossible.  Notice  that  the  profiles  given  in  Fig.  34 
are  the  exact  solutions  of  the  PFT  equations  (45-46)  associated  with  the  thermal 
equilibrium,  therefore,  they  provide  the  exact  relation  for  the  chemical  potential 
jLi  =  constONQX  all  the  domain  of  solution.  Within  the  three-phase  interface,  this 
relation  holds  as  well,  while  the  concentration  and  phase  distributions  are  clearly 
not  uniform.  Thus,  it  is  natural  to  assume  that  the  changes  observed  in  the  phase 
diagram  obtained  by  solving  PFT  equations  are  due  to  the  account  of  interfacial 
energy. 

Thus,  we  see  that  the  standard  PFT  model  (45-47)  formulated  above 
reproduces  the  general  properties  of  phenomena  associated  with  the  formation  of 
DZ  and  CM  near  the  eutectic  point,  in  particular,  in  the  systems  with  amorphous 
Si.  As  for  the  experiments  with  crystalline  Si  films,  they  may  reveal  more 
complicated  properties  of  the  DZ  [20].  Apparently,  to  describe  these  latter  more 
adequately,  the  detailed  description  of  component  structures  occuring  in  DZ  is 
needed,  eventually,  by  using  the  multi-phase  version  of  PFT  [28,  29].  At  the  same 
time,  the  basis  for  the  consideration  is  the  assumption  is  that  the  relatively  wide  DZ 
is  associated  with  the  equilibrium  state  of  complete  decomposition  with  high 
interfacial  solid-solid  energy,  and  the  crucial  for  the  model  parameter,  the  DZ 
width,  still  remains  of  the  same  order. 
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The  results  of  simulations  suggest  that  the  three-phase  states  and  the 
assoeiated  properties  of  euteetie  phase  diagram  would  be  rather  eommon  for 
euteetie  systems.  Another  phenomenon,  whieh  ean  be  viewed  as  the  manifestation 
of  the  formation  of  three-phase-states,  is  the  miero-separation  of  phases  observed 
within  the  liquid  euteetie  above  the  euteetie  point  [13,14].  The  experiments  elearly 
demonstrate  the  strong  deviation  from  the  regular  liquid  strueture  in  liquid  euteetie 
even  on  maeroseopie  seales,  especially,  after  the  long-term  (during  several  hours) 
keeping  in  isothermal  conditions. 

Conclusions 

By  means  of  computer  simulations  based  on  the  standard  PFT  model  with 
the  miscibility  gap  in  solid  state,  the  following  issues  were  examined:  (1)  structure 
formation  in  binary  LaBg -ZrBj  eutectic  as  crystallized  from  the  overcooled  melt; 
(2)  formation  of  ordered  fiber  structures  in  directional  crystallization  of 
LaBg-ZrB2;  (3)  interfacial  phenomena  in  binary  eutectic:  the  formation  of 
diffusion  zone  and  contact  melting  in  Cu-Si  system. 

The  simulations  as  a  whole  reproduce  the  basic  features  of  eutectic 
structures  observed  in  experiments  with  boride-boride  ceramics,  -  the  spatial 
component  segregation  and  formation  of  structure  from  glass,  the  complete  system 
decomposition  in  the  process  of  relaxation  to  the  thermodynamical  equilibrium, 
formation  of  colonies  in  directional  crystallization.  Further  progress  in  description 
of  structure  formation  in  boride  composites  should  be  expected  with  refining  the 
experimental  parameters  needed  for  formulation  of  PFT  models,  and  with 
implementation  of  more  detailed  3D  models  [28,29],  which  take  into  account  the 
differences  in  the  crystal  structure  of  the  components. 

The  simulations  of  the  interfacial  phenomena  in  Cu-Si  predict  the  existence 
of  equilibrium  three-phase  (solid-liquid-solid)  states  above  the  eutectic  temperature 
describing  the  complete  decomposition  with  a  liquid  interface.  The  results  of 
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simulations  demonstrate  the  role  of  the  interfaeial  energy  as  the  possible 
meehanism  for  the  phenomena  of  eontaet  melting  and  formation  of  diffusion  zone 
observed  in  the  experiments  with  binary  metal-silicon  systems.  In  addition,  the 
results  obtained  provide  the  explanation  of  the  phenomenon  of  phase  separation  in 
liquid  eutectic. 
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II.  MICROSCALE 


5.  Calculation  of  physical  and  mechanical  characteristics  of  MeB2  (Me  -  Ti, 
Zr,  Hf)  BORIDES;  LaBg  AND  composites  in  the  LaBg  -  MeB2  systems  on  the 
basis  of  a  priori  pseudopotential 

The  promising  prospects  for  the  use  of  materials  on  the  basis  of  LaBg  -  MeB2 
in  various  branches  of  engineering  cause  high  interest  in  determining  their  strength 
characteristics  for  which  data  are  virtually  absent  [30-32]  especially  for  composite 
materials.  Therefore,  the  issue  of  calculating  the  strength  characteristics  of 
composites  LaB6  -  MeB2  (Me  -  Ti,  Zr,  Hf)  with  eutectic  composition  is  high  on 
the  agenda. 

In  the  previous  project  (HOARD  N  P-273)  a  model  was  developed  [33], 
which  allowed  to  construct  thermodynamic  potentials  of  quasi-binary  eutectic 
systems  LaB6-MeB2,  and  calculate  for  them  the  concentrations  of  the  components 
and  melting  point  temperature  in  the  eutectic  point. 

The  energy  of  the  electron-ion  system  was  calculated  using  method  of  a 
priori  pseudopotential  [34]  and  concept  of  intermolecular  interaction  [35]. 
Intermolecular  interactions  by  nature  are  close  to  interatomic  interactions  and  are 
described  by  the  same  kind  of  potentials  as  the  interatomic.  The  concept  of 
intermolecular  interaction  becomes  useful  when  it  comes  to  interaction  of  the 
components  in  the  alloy,  especially  when  the  components  are  made  up  of  different 
atoms  (LaB6  or  MeB2). 

In  the  method  of  pseudopotentials  the  energy  of  electron-ion  system 
attributed  to  the  unit  cell  (molecule)  is  the  sum  of  the  energy  of  electrons  and  ions 
included  in  its  structure.  Part  of  this  energy,  which  depends  on  the  type  of  crystal 
lattice  energy  is  responsible  for  the  pairwise  intermolecular  interaction.  In 
calculating  the  energy  of  interaction  between  heterogeneous  molecules  the  virtual 
crystal  model  approximation  is  accepted.  For  this  purpose  Helmholtz 
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thermodynamic  potential  is  built,  which  depends  on  the  temperature  and 
concentration  of  the  components.  From  the  extreme  condition  of  the 
thermodynamic  potential,  were  calculated  characteristic  parameters  of  the  eutectic 
(the  concentration  of  the  component  and  the  melting  point  temperature)  of  LaB^  - 
MeB2  system,  values  of  which  are  in  good  agreement  with  experiment  [33,  36]. 

Further  progress  of  this  work  consists  in  investigation  of  physic-mechanical 
properties  of  LaBg  -  MeB2  systems  under  extreme  conditions  (high  temperature 
and  deformation).  Taking  into  account  complexity  of  the  formulated  problem,  we 
approach  its  solution  on  stage-by-stage  basis: 

1.  Calculation  of  mechanical  characteristics  of  components  constituting  the 
composite.  For  calculation  of  these  characteristics  for  investigated  LaB6-MeB2 
systems  under  uniaxial  deformations,  we  use  the  rule  of  mixtures  [37-40]  as  a  first 
approximation. 

2.  The  problem  of  calculation  of  thermal  expansion  coefficients  of 
composites  with  eutectic  composition  becomes  even  more  acute  when  these 
materials  are  considered  for  use  under  extreme  operational  conditions.  Difference 
in  thermal  expansion  coefficients  of  components  leads  to  emergence  of  internal 
stresses,  neglect  of  which  makes  control  of  mechanical  properties  less  effective. 
This  reduces  range  of  application  of  such  materials.  That  is  why  preliminary 
estimate  of  such  significant  material  characteristic  as  thermal  expansion 
coefficient,  especially  in  absence  of  both  experimental  and  theoretical  data  for  this 
parameter  [41],  is  expedient. 

3.  To  find  dependence  of  mechanical  characteristics  (deformation,  elasticity 
modulus,  theoretical  strength)  from  temperature,  one  has  to  know  energy  of 
electron-ion  system  of  materials  at  different  temperatures.  In  pseudopotential 
method  this  means  necessity  to  determine  volume  change  of  elementary  cell  of 
borides  at  temperatures  different  from  zero,  that  is  obtain  explicit  dependence  of 
full  energy  of  the  system  from  lattice  parameters  in  a  given  interval  of 
temperatures.  For  this,  we  use  a  physical  model  developed  to  calculate  thermal 


Distribution  A:  Approved  for  pubiic  reiease;  distribution  is  uniimited. 


expansion  coefficient,  and  obtain  the  dependence  of  lattice  parameters  (volume) 
from  temperature  [42]. 

4.  On  the  basis  of  a  computational  experiment,  a  model  of  the  formation  of 
interphase  boundaries  in  the  quasi-binary  eutectic  systems  is  proposed.  The  values 
of  the  energy  of  the  interface  between  the  components  in  the  boride  and  metal- 
ceramic  quasi-binary  composites  with  eutectic  composition  were  calculated.  The 
values  of  the  energy  of  the  interface  at  different  temperatures  for  LaB6-TiB2,  LaB^- 
ZrB2  systems  [43]  were  estimated. 

5.  Calculation  of  the  mechanical  characteristics  of  quasi-binary  eutectic 
composites  with  account  of  inter-component  interaction  at  the  interface  between 
components  [44]. 

Calculation  of  theoretical  strength  and  Young's  modulus  ofMeB2  (Me  -  Ti,  Zr) 
borides,  LaB^  composites  with  eutectic  composition  in  LaB^-MeB2  systems 

To  calculate  mechanical  characteristics  of  the  material  on  the  micro  level  the 
value  of  the  total  energy  per  atom  (molecule)  of  the  material  is  used  [33]. 

Energy  (at  constant  volume  of  the  molecule)  is  a  sum  of  the  electrostatic 
energy  and  the  energy  of  band  structure,  which  can  be  represented  as  a  sum  of 
paired  intermolecular  potentials.  The  change  in  energy  along  the  lattice  parameter 
(or  its  projection)  parallel  to  the  axis  of  deformation  is  directly  proportional  to  the 
theoretical  strength  of  the  material.  Under  uniaxial  deformation  along  the  z-axis 
strength  is  determined  through  the  energy  of  the  electron-ion  system  (U)  per  unit 
area  of  atomic  plane  in  the  unit  cell,  located  perpendicular  to  the  axis  of 
deformation  [45] 

(7^  =(d  UldQ^  )/(S-  d),  (49) 

where  S  -  surface  of  atomic  plane  in  atomic  cell,  -  relative  deformation,  d  - 
lattice  parameter  or  its  projection  parallel  to  deformation  axis. 

At  very  small  strains  stress-strain  curve  is  close  to  linear  law  and  then: 
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G  -  E-Adido. 


(50) 


In  the  case  of  composites  for  calculation  of  Young's  modulus  it  is  necessary 
to  make  some  simplifications.  Assume  that  the  reinforcing  phase  (fibers  or  plates) 
is  distributed  homogeneously  and  its  basic  structural  elements  are  parallel  to  each 
other  throughout  the  matrix.  In  most  cases,  elastic  modulus  of  the  composite 
material  can  be  calculated  by  the  rule  of  mixtures.  However,  to  use  this  it  is 
necessary  to  clarify  the  type  of  loading.  Consider  the  homogeneous  strain  state 
when  the  reinforcer  is  located  along  the  direction  of  the  load.  Then,  for  the  first 
stage  (region)  of  deformation,  where  the  matrix  and  the  reinforcer  deform 
elastically,  for  the  theoretical  strength  of  the  composite  it  can  be  written  [46] 

Gr  —  G^  +  dg  Gg  —  d^  +  EgSgdg  (51) 

where  Ga  -  stress  in  the  matrix  at  fiber  fracture  strain,  gb  -  strength  limit  of  fibers 
under  tension;  d^,dg  -  volume  fraction  of  components  (Table  6),  e^,ej^  - 
deformations. 

Table  6.  Volume  fraction  of  components  in  composite 


CncxeMa 

4 

LaB6  -  TiB2 

0.893 

0.107 

LaB^  -  ZrB2 

0.837 

0.163 

LaBa  -  HfB2 

0.873 

0.126 

The  dependence  of  stress  from  strain  for  components  obtained  by  the 
formula  (49)  are  shown  in  Figure  39.  Less  ductile  are  MeB2,  maximum 
deformation  for  these  diborides  is  about  0.1,  while  for  LaBe-O.lh.  TiB2  has  greater 
strength  compared  with  ZrB2  and  HfB2  diborides,  strength  values  of  which  are 
almost  identic. 
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The  strength  of  the  euteetie  composite  with  brittle  hardening  phase,  oriented 
along  the  axis  of  tension  (Fig.  40)  can  be  obtained  by  substituting  in  (51) 
£^  =  £g=  £^^^ ,  where  Smax  -  maximum  strain  of  fragile  phase  in  tension.  It  is 
believed  that  the  fracture  of  the  composition  occurs  as  a  result  of  simultaneous 
fracture  of  equally  strong  fibers  in  one  section  when  the  maximum  deformation 
£max  is  rcachcd.  As  TiB2,  ZrB2  and  HfB2  deform  almost  elastically,  then  the  use  of 
mixture  rule  for  system  LaB6-MeB2  when  calculating  mechanical  characteristics  is 
an  acceptable  approximation. 

The  results  of  computational  experiment  to  determine  elastic  moduli  of  the 
components  and  composite  for  investigated  systems  by  the  formulas  (49-51)  are 
presented  in  tabular  form. 


Table  7.  The  elastic  moduli  of  the  components  and  composite  materials  in  the 
system  LaB6-MeB2  and  their  experimental  value  in  GPa 


Phase 

E  (calc.) 

E  (exp.) 

LaBg 

505.27 

478.73  [47]  320" 

TiB2 

608.84 

540.53  [47]  545" 

HfB2 

523.97 

479.71  [47] 

ZrB2 

534.67 

495.80  [47]  430" 

LaB6-TiB2 

516.11 

- 

LaB6-HfB2 

507.63 

- 

LaB6-ZrB2 

508.46 

430-450" 

(^Private  communication  A.B.  Lyashchenko) 


The  obtained  values  of  Young's  modulus  are  close  to  their  experimental 
values  within  acceptable  limits  (maximum  relative  error  of  ~  7%),  taking  into 
account  that  in  the  calculations  we  considered  only  perfect  single  crystals,  while 
real  crystals  always  have  lower  values. 
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a,  GPa 


Fig.  39.  Dependence  of  strength  from  strain  for  borides  at  T  =0  K 
a,  GPa 


Fig.40.  Dependence  of  theoretical  strength  (at  T  =  0  K)  1  -  TiB2 , 2  -  LaB6  and 
system  3  -  LaB6  -  TiB2  with  eutectic  composition,  calculated  by  mixture  rule 
taking  into  account  their  volumetric  fraction  from  deformation. 
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Young's  modulus  and  theoretical  strength  (Table  7)  of  borides  of  transition 
metals  decreases  with  increasing  atomic  number  of  metals  in  the  group  similar  to 
the  case  of  experimental  measurements  of  microhardness. 

Note  that  the  equations  of  "rule  of  mixtures "  is  adopted  in  the  theory  of 
artificial  composites  with  fiber  reinforcement  (Kurz  W.,  Sahm  PR  Gerichtet 
erstarrte  eutektische  werkstoffe  -  Springer-Verlag,  Berlin  / Heidelberg.-  1975.  -p.- 
272).  It  is  assumed  that  the  interaction  between  the  phases  is  absent.  This 
assumption  is  valid  since  in  artificial  composite  sizes  of  particles  of  phase 
components  are  large,  and  bonding  at  interfaces  is  weak.  On  the  contrary,  eutectic 
phases  exert  significant  influence  on  each  other  at  the  interface.  Thus,  the  rule  of 
mixtures  in  application  to  eutectic  composites  can  be  viewed  as  a  first 
approximation.  Hence,  further,  mathematical  models  are  developed  that  can  be 
used  to  calculate  mechanical  properties  of  quasi-binary  composites  taking  into 
account  interaction  between  phases  at  the  interface. 

Calculation  of  linear  thermal  expansion  coefficients  of  MeB2  (Me  -  Ti,  Zr) 
borides,  LaB^  and  eutectic  composites  in  LaBg-MeB2  systems 

The  scatter  of  experimental  data  on  the  coefficients  of  thermal  expansion  for 
borides  of  transition  metals  for  the  same  compounds  sometimes  can  be  as  high  as 
70-100%  [47].  This  is  due  to  he  lack  of  sufficiently  pure  single  crystals.  The 
presence  of  a  small  percentage  of  impurities  can  greatly  alter  physical  properties  of 
the  studied  systems.  In  addition,  while  for  borides  (MeB2,  LaBg)  thermal  expansion 
coefficients  are  experimentally  determined,  and  this  variation  is  established,  then 
for  LaBa  -  MeB2  system  published  experimental  data  are  not  available.  Therefore, 
calculation  of  thermal  expansion  coefficients  for  borides  (MeB2,  LaBa),  as  well  as 
quasi-binary  systems  LaB6-MeB2  eutectic  composition,  based  on  first  principles,  is 
a  matter  of  topical  interest. 
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It  should  be  noted  that  the  method  of  preparing  the  samples,  the  phase 
eomposition,  the  presenee  of  impurities  and  the  measurement  method  have 
significant  influence  on  the  values  of  the  physical  characteristics  of  borides. 
Theoretical  calculations  consider  ideal  single  crystals  and  it  usually  gives 
somewhat  higher  values  of  the  calculated  parameters  compared  with  the 
experiment. 

One  of  possible  ways  of  calculating  the  coefficient  of  thermal  expansion  is  a 
method  based  on  determination  of  the  dependence  of  crystal  total  energy  on  the 
lattice  parameters  at  different  temperatures. 

The  total  energy  of  the  crystalline  material  can  be  represented  as  the  sum  of 
the  energies  of  the  electron-ion  system  at  T  =  0,  (T  -  absolute  temperature)  and  the 
energy  of  the  thermal  vibrations  of  the  ions  at  temperature  T  ^  0.  To  calculate  the 
energy  of  the  electron-ion  system  at  T  =  0  we  use  the  method  of  a  priori 
pseudopotential  [34].  The  energy  of  the  thermal  vibrations  of  the  investigated 
materials,  which  have  complex  structure,  we  calculate  within  the  framework  of 
Einstein's  model  [48],  which  gives  reliable  results  in  contrast  to  Debye  model  (that 
works  quite  well  for  simple  crystal  lattices  [49]  and  leads  to  unstable 
computational  schemes  for  complex  structures). 

The  total  energy  of  the  electron-ion  system  at  T  =  0  within  the  framework  of 
the  pseudopotential  method  in  the  second-order  of  perturbation  theory,  can  be 
written  as  [33] 


U  =  u,+u,,  (52) 

where  [/q  depends  on  the  volume  of  the  unit  cell.  Energy  Uq  includes  kinetic, 

exchange-correlation  energy  and  energy  of  a  free  electron  gas,  as  well  as 
corrections  to  the  energy  in  the  first  order  of  perturbation  theory  with  respect  to 
pseudopotential.  The  second  term  (second-order  correction)  represents  energy  of 
band  structure  and  electrostatic  energy,  and  can  be  represented  as  a  sum  of  pair 
interatomic  potentials 
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(53) 


ly  i,  j 


where  R.  -Rj  -  distance  between  atoms  i  and  j,  N  -  number  of  atoms  in  a 

representative  volume  [50].  The  method  of  calculating  the  energy  of  the  electron- 
ion  system  in  the  framework  of  the  pseudopotential  method  is  given  in  [33]. 

According  to  Einstein's  model  [48],  atoms  in  the  crystal  lattice  occillate  with 
the  same  frequency,  which  is  proportional  to  material  stiffness.  The  average  energy 
of  the  lattice  vibrations  is  defined  by  [48] 


q 


flQ) 

_ q _ 

exp(/i<y^  /  kT)  - 1  ’ 


(54) 


where  summation  is  done  over  all  types  of  vibrations,  and  co^  -  oscillation 
frequency  (according  to  Einstein  model  (o^=co  for  all  q).  In  the  case  of  the 

structure  with  n  atoms  in  a  unit  cell,  there  are  3iVn  oscillations  {N  -  number  of  cells 
in  the  material  volume)  and  the  oscillation  frequency  is  determined  by  the  relation 

CO  =  ^lla  I M  .  (55) 


Here  a*  -  force  constant,  which  is  determined  by  the  second  derivative  of  the 
energy  of  interatomic  interaction  with  respect  to  space  variable 


* 


a 


( ;)2 


(56) 


Pair  interatomic  potentials  (53)  for  the  investigated  materials  were  constructed  in 
[33]. 

In  calculating  the  energy  of  the  electron-ion  system  in  the  second-order  of 
perturbation  theory  with  respect  to  the  pseudopotential,  we  use  harmonic 
approximation.  However  the  application  of  this  approximation  in  the  lattice 
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dynamics  is  not  sufficient  for  calculation  of  the  coefficient  of  thermal  expansion 
from  the  “first  principles”.  At  the  same  time,  we  know  that  the  calculation  of  the 
vibrational  part  of  corresponding  thermodynamic  functions  taking  into  account 
anharmonicity  is  hopelessly  complex  [51].  This  difficulty  can  be  avoided  if  we 
make  simplifying  assumptions  about  the  dependence  of  the  oscillation  frequencies 
on  temperature  and  volume  using  the  “quasi-harmonic”  approximation  [51].  In  this 
approach,  we  first  calculate  force  constants  at  T  =  0.  The  energy  of  thermal 
vibrations  (54),  calculated  at  low  temperatures  (rj=Ar;  AT  ^0)  and  force 

constants  at  T  =  0  is  added  to  the  value  of  the  energy  U.  At  the  subsequent 
steps  in  the  calculation  of  the  total  energy  within  the  harmonic  approximation,  we 
use  new  power  constants,  and,  therefore,  frequencies  dependent  on  volume. 

Schematically,  the  solution  of  the  formulated  problem  can  be  presented  as 
follows: 

-  From  the  minimum  of  the  energy  of  the  electron-ion  system  U(a,c)  lattice 
parameters  (ao,  Cq)  and  force  constants  in  the  equilibrium  state  of  the  crystal  are 
determined. 

-  Full  energy  U*  =1/^+1/^  of  the  system  at  the  temperature  T  ^0  is  calculated. 

-  From  the  equation  U*^^(aQ,CQ,T)  =  U new  values  of  lattice  parameters  are 
found. 

-  From  formula  (56),  at  the  obtained  values  of  lattice  parameters,  the  new  values 
of  the  force  constants  are  calculated. 

As  a  result,  we  obtain  the  dependence  of  lattice  parameters  from 
temperature. 

Coefficient  of  linear  thermal  expansion  a.  is  determined  from  the  relation 
[30] 


A  a,. 

O',  = — ^ 
a. 


T 


(57) 


{tti  -  lattice  parameter,  A  a,  -  change  of  the  lattice  parameter  at  the  temperature  7). 
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For  LaB^  with  cubic  structure  lattice  parameter  is  ai=  a,  and  for  hexagonal 
MeB2  ai=  {a\  c). 

Fig.  41-43  shows  the  temperature  dependence  of  relative  elongation,  and 
Tables  8  and  9  -  corresponding  values  of  the  coefficients  of  thermal  expansion  for 
the  materials  studied. 

Values  of  ««  are  almost  constant  in  the  whole  range  of  temperatures  studied 
for  ZrB2,  and  in  the  case  of  TiB2  —  monotonously  increase  up  to  temperature  of 
1000  K  and  then  stabilize.  Coefficient  of  thermal  expansion  of  LaB6  (Table  6)  has 
a  jump  in  the  temperature  range  from  750  K  to  1000  K.  This  interesting  result  has 
quite  simple  physical  sense.  The  basis  for  the  calculation  of  the  physical 
characteristics  of  crystals  from  first  principles  is  the  type  of  crystal  potential.  This 
potential  is  used  for  construction  of  unit  cell  of  the  crystal  under  study  from 
corresponding  atoms,  cell  type  is  determined  based  on  the  minimum  energy  of 
electron-ion  system,  as  well  as  the  parameters  of  the  crystal  lattice.  According  to 
the  results  of  calculations,  the  energy  of  interaction  of  boron  atoms  (complex  Bg) 
in  the  structure  of  LaBe  has  a  deeper  minimum  than  the  energy  of  interaction 
between  atoms  of  La  and  B,  or  atoms  of  La-La.  Therefore,  in  the  initial  stage  of 
raising  the  temperature  to  750  K  it  is  La  atoms  that  are  displaced  from  the 
equilibrium  position.  At  higher  temperatures  Be  volume  complex  begins  to 
participate  in  the  expansion,  which  leads  to  a  jump  in  the  coefficient  of  thermal 
expansion. 


Distribution  A:  Approved  for  pubiic  reiease;  distribution  is  uniimited. 


Table  8.  Coefficients  of  thermal  expansion  ««  (in  K'^  units)  for  TiB2,  ZrB2,  LaBg 
and  their  averaged  experimental  values  [47]  (in  brackets) 


T,K 

aa  (TiB2) 

aa  (ZrB2) 

ac  (LaBe) 

(4.5). 

(5.9). 

(6.4±0.5). 

300 

5.503 

6.144 

7.597 

500 

5.540 

6.144 

7.597 

750 

5.590 

6.147 

7.597 

1000 

5.70 

6.145 

7.850 

1500 

5.655 

6.147 

7.851 

2000 

5.701 

6.149 

7.850 

2500 

5.699 

6.146 

7.850 

2750 

5.762 

6.145 

7.801 

In  contrast  to  the  cubic  LaBe,  transition  metal  diborides  have  chracteristic 
hexagonal  structure,  so  the  task  is  set  of  calculation  of  the  coefficient  of  thermal 
expansion  also  in  the  direction  perpendicular  to  the  basal  plane  (direction  c).  To 
calculate  the  coefficient  of  thermal  expansion  in  the  direction  of  c,  we  propose  to 
use  the  energy  of  interaction  between  the  basal  planes,  whose  value  for  the  systems 
studied  is  calculated  in  [38].  At  high  temperatures,  the  distance  between  the  basal 
planes  increases,  which  leads  to  a  change  in  the  interaction  energy  between  the 
planes. 
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Table  9  Coefficients  of  thermal  expansion  for  TiB2  and  ZrB2,  in  brackets  their 
averaged  experimental  values  are  given  [47]. 


T,K 

ac-lO^K^  (TiB2) 

a.-lO^r^  (ZrB2) 

(8.65-11.20) 

(6.78-7.65) 

300 

10.231 

7.589 

500 

10.276 

7.589 

750 

10.472 

7.659 

1000 

10.706 

7.750 

1500 

10.840 

7.808 

2000 

11.072 

7.790 

2500 

10.894 

7.684 

2750 

10.267 

7.313 

3000 

10.012 

7.083 

Typically,  the  components  of  the  composition  (the  matrix  and  the  reinforcer) 
have  different  coefficients  of  thermal  expansion.  Approximately  linear  expansion 
coefficient  of  the  composite  is  determined  by  the  relation  [52] 


where  and  E’g  are  Young  moduli  of  components,  and  E  -  average  modulus  for 
composition 


E  =  +  S,E,  .  (59) 

Young's  moduli  are  determined  through  the  stress  tensor  [53]. 

The  coefficient  of  thermal  expansion  of  the  composite  in  a  direction  a  is 
calculated  by  the  formula  [52] 

~  ^  A  ‘  ^ A  ^ B  ‘  ^ B 
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Given  that  the  composites  have  a  fiber  structure  and  fiber  diameter  is  much  less 
than  its  length,  formula  (60)  was  applied  [52].  On  the  other  hand  thermal 
coefficient  for  direction  c  was  also  calculated  by  the  (60).  The  results  showed  that 
the  coefficients  calculated  by  the  formulas  (56)  and  (60)  do  not  differ  much  from 
each  other. 

The  calculation  results  for  the  composites  are  presented  graphically: 

1)  relative  elongation  of  the  composite  with  eutectic  composition  in  the 
direction  perpendicular  to  the  solidification  front  (Figure  41)  (direction  c  of 
the  hexagonal  lattice); 

2)  dependence  of  the  thermal  expansion  coefficient  from  temperature  (Fig.  42). 

Results  of  computational  experiment  on  the  Ab-initio  calculation  of  the 
coefficient  of  thermal  expansion  of  boride-boride  eutectic  composites  show  [41]: 
Investigated  composites  have  pronounced  anisotropy  of  properties,  thermal 
expansion  coefficients  and  are  very  different  from  each  other. 

In  the  temperature  range  750  K<T  <\  100  K  thermal  expansion  coefficients 
and  ttc  for  the  studied  composites  LaB6-MeB2  change  abruptly,  due  to 
abrupt  change  in  CTE  of  LaBg. 

The  average  values  of  the  thermal  expansion  coefficients  ««  for  composites 
and  LaB6-TiB2  are  almost  identical,  and  ac  for  LaB6-ZrB2  over  the  entire 
range  of  temperature  is  higher  than  that  for  LaB6-TiB2. 

The  coefficients  of  thermal  expansion  (««  and  a^)  for  TiB2  and  ZrB2,  as  well 
as  for  LaBg  are  calculated  with  the  help  of  the  proposed  model  using  a  priori 
pseudopotential  method,  and  are  within  the  scatter  of  experimental  data. 
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Aa/a  *  10 


LaB, 


Fig.  41.  Temperature  dependence  of  the  relative  elongation  along  a  (for  MeB2  a 
parameter  of  the  basal  plane) 


Ac/c  *  1(f 


TiB^ 


Fig.42.  Temperature  dependence  of  the  relative  elongation  in  the  direction 
perpendicular  to  the  basal  planes  for  components  with  a  hexagonal  lattice. 
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Fig.  43.  Dependence  of  the  thermal  expansion  coefficients  and  «c  on  the 
temperature  of  composites  with  eutectic  composition. 

Since  LaB^  has  a  structure  of  CaBg  type  (simple  cube  of  metal  atoms 
centered  by  octahedron  of  boron  atoms),  its  thermal  properties  in  the  direction  a  or 
c  are  identical.  The  anisotropy  of  the  thermal  coefficient  of  the  composite  is 
connected  with  the  second  components  MeB2.  Diborides  of  transition  metals  have 
a  hexagonal  structure  of  AIB2  type,  in  which  boron  atoms  form  graphite-like  grids 
perpendicular  to  the  axis  Z,  and  the  entire  structure  -  is  a  consequent  alternation  of 
hexagonal  layers  of  metal  atoms  arranged  in  a  hexagonal  close-packed  lattice  sites 
with  a  small  da  ratio  and  layers  of  boron  atoms  forming  a  hexagonal  two- 
dimensional  grid.  According  to  calculations,  lattice  parameters  of  TiB2  and  ZrB2 
are  a  =  3.05  A;  c  =  3.24  A  and  a  =  3.12  A;  c  =  3.49  A  [33].  On  the  basal  plane 
lattice  parameters  for  two  diborides  are  close  in  value,  which  can  not  be  said  about 
the  parameters  in  the  direction  perpendicular  to  the  basal  plane.  Physical- 
mechanical  properties  of  materials,  including  thermal,  are  based  on  the  interatomic 
interactions,  and  the  latter  depend  on  the  interatomic  distances.  This  fact  is 
reflected  in  the  results  of  the  calculations,  which  supports  the  adequacy  of  the 
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models  chosen.  The  composition  of  LaB6-ZrB2  and  LaB6-TiB2  composites  contain 
70%  and  75%  of  LaB6,  respectively,  and,  hence,  their  characteristics  calculated  by 
the  formulas  (58)  mainly  reflect  the  thermal  properties  of  this  component. 
Therefore,  in  the  same  temperature  range  from  750  K  to  1000  K,  the  abrupt 
increase  in  the  coefficient  of  thermal  expansion  of  these  composites  is  observed. 

The  values  of  the  coefficients  of  thermal  expansion  of  LaB6-ZrB2  and  LaBe- 
TiB2  composites  obtained  using  formulas  (58  -  60)  are  less  accurate  compared  with 
similar  calculations  for  the  components,  as  in  the  case  of  composites  at  this  stage 
the  influence  of  the  contact  borders  of  the  two  components  is  not  taken  into 
account. 


Temperature  dependence  of  the  theoretical  strength  of  borides  and  quasibinary 
boride  eutectic  systems 

With  increasing  temperature,  the  level  of  mechanical  properties  that 
characterize  strength  of  the  material  decreases  and  the  level  of  plastic  properties  — 
increases.  Many  of  the  materials  at  elevated  temperatures  undergo  physico¬ 
chemical  transformations  (precipitation  of  hardening  phases,  coagulation  and 
dissolution  of  phases,  oxidation,  especially  of  grain  boundaries,  etc.).  Depending 
on  the  nature  of  these  transitions  noticeable  deviations  in  the  temperature  curves  of 
strength  and  ductility  can  be  observed,  as  well  as  change  in  the  character  of 
fracture. 

Physical  properties  that  determine  material  characteristics  at  high 
temperatures  are  energy  of  interatomic  bonding  and  melting  point  temperature  of 
compounds.  Therefore,  the  alloys  based  on  refractory  metals  weaken  less  with 
temperature  than  the  alloys  having  low  melting  point  and  energy  of  bonding  forces. 

One  of  the  main  requirements  for  the  composites,  the  end  product  made 
from  which  are  used  in  extreme  conditions  —  in  internal  combustion  engines,  gas 
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pipes,  aircraft  jet  engines  —  it  is  their  performance  at  high  temperatures.  This  is 
what  explains  the  interest  to  the  study  of  variation  in  the  mechanical  properties  of 
composites  as  a  function  of  operational  temperature.  Effect  of  high  temperature  on 
the  mechanical  properties  of  eutectic  boride  and  ceramic  composites  and  functional 
materials  is  presented  in  the  work  by  S.S.  Ordanyan  et  al.  [30].  It  investigated  the 
temperature  dependence  of  the  bending  strength  for  metal-ceramic  and  boride 
systems  and  showed  that  in  the  studied  systems  of  eutectic  composition  high 
strength  is  retained  at  temperatures  up  to  0.8  times  the  eutectic  temperature  (Tev). 

Dependence  of  the  mechanical  properties  (microhardness,  flexural  strength) 
on  the  temperature  in  the  range  from  25  to  1600  °C  have  been  studied  in  [31,  32, 
54]  for  LaB6,  TiB2  and  LaB6-TiB2  system,  where  it  is  assumed  that  the  flexural 
strength  in  the  directionally  reinforced  LaB6-TiB2  composite  at  high  temperatures, 
mainly  depends  on  the  ductility  of  the  LaB6  matrix  and  TiB2,  which  is  a  rather 
obvious  fact.  This  is  a  consequence  of  the  formula  that  determines  the  strength  of 
the  system  by  the  rule  of  mixtures,  where  the  decisive  role  is  played  by  the 
maximum  strain  of  diborides. 

From  the  works  on  computational  experiments  related  to  the  study  of 
borides,  one  may  point  out  [55],  where  with  the  use  of  density  functional  theory, 
TiB2  bulk  modulus  was  calculated  and  the  compount  was  classified  as  a  superhard 
material  (bulk  modulus  of  about  260  GPa).  In  [56],  on  the  basis  of  the  Hartree- 
Fock  method  and  density  functional  theory,  anisotropic  elastic  constants  for  TiB2 
were  calculated.  However,  these  studies  did  not  consider  the  dependence  of 
mechanical  properties  on  the  temperature  for  transition  metal  borides,  and 
especially  for  FaB6-MeB2  systems.  This  is  due  to  the  complexity  of  the  calculation 
of  thermal  energy  of  the  system.  With  increasing  temperature,  the  thermal 
vibrations  of  the  atoms  are  not  subject  to  the  harmonic  approximation.  Accounting 
for  non-harmonic  members  is  a  rather  difficult  task  in  the  framework  of  quantum- 
mechanical  calculations,  still  a  question  of  dependency  of  the  mechanical 
characteristics  on  the  temperature  remains  relevant. 
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In  [41],  where  the  calculation  of  the  linear  coefficient  of  thermal  expansion 
(LCTE)  of  borides  from  ab  initio  methods  was  done,  the  dependence  of  LCTE  on 
the  temperature  was  revealed,  that  is  a  purely  quantum-mechanical  effect.  The 
model  developed  in  [41]  using  quasi-harmonic  approximation  makes  it  possible 
within  the  framework  of  quantum-mechanical  calculations  to  study  the  influence  of 
temperature  on  the  tensile  strength  of  borides  and  boride-boride  eutectic  systems. 

To  identify  the  dependence  of  mechanical  characteristics  on  temperature,  it 
is  necessary  to  be  able  to  calculate  the  energy  of  the  electron-ion  system  of  the 
materials  at  different  temperatures.  Within  the  framework  of  the  pseudopotential 
method  it  means  to  find  the  change  of  volume  of  the  unit  cells  of  borides  at 
temperatures  different  from  zero,  i.e.,  obtain  the  explicit  dependence  of  the  total 
energy  from  the  lattice  parameters  or  volume  at  nonzero  temperature. 

To  solve  this  problem  it  is  necessary: 

-  Erom  the  minimum  energy  of  the  electron-ion  system  U(Q.)  determine  the 
volume  of  the  lattice  (Q)  in  the  equilibrium  state  of  the  crystal. 

-  Calculate  the  total  energy  of  the  system  U*  =Uq+Uj  at  nonzero 
temperatures  (  r  ^  0). 

-  Erom  equation  =  U(Qi)  determine  the  new  values  of  the  volume 

of  the  crystal  lattice. 

Einally,  we  obtain  dependence  of  the  volume  from  the  temperature. 

The  procedure  for  the  calculation  of  the  energy  of  the  electron-ion  system  in 
the  framework  of  the  pseudopotential  method  is  given  in  [33,  56  -  58].  The  results 
of  computational  experiment  to  determine  the  dependence  of  the  volume  of  the 
unit  cells  from  the  temperature  are  presented  in  tabular  form  (Table  10). 
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Table  1.5.  Calculated  values  of  the  volume  of  the  unit  cells  of  investigated  borides 
depending  on  the  temperature 


T,  K 

Q  (TiB2 ) 

Q  (ZrB2) 

Q  (LaBe ) 

0 

173.238 

207.810 

483.770 

300 

174.108 

209.054 

487.219 

500 

174.696 

209.782 

489.541 

750 

175.432 

210.692 

492.444 

1000 

176.168 

211.602 

495.346 

1500 

177.640 

213.421 

501.151 

2000 

179.113 

215.241 

506.956 

2500 

180.585 

217.061 

512.761 

2750 

181.347 

217.987 

515.663 

Under  uniaxial  deformation  along  the  z  axis  theoretical  strength  is 
determined  by  the  energy  of  the  electron-ion  system  U  per  unit  area  of  atomic 
plane  in  the  unit  cell  perpendicular  to  the  axis  of  deformation 


du^ 

C  dCy 


(61) 


where  -  relative  deformation,  c  -  lattice  parameter  or  its  projection  parallel  to  the 
axis  of  deformation.  For  each  value  of  the  volume  (corresponding  to  a  given 
temperature)  the  energy  of  the  electron-ion  system  is  calculated.  Theoretical 
strength  is  calculated  with  formula  (61)  by  standard  methods  [39].  As  a  result,  the 
computational  experiment  found  that  at  temperatures  close  to  2000  K  the  value  of 
the  absolute  deformation  for  diborides  is  increased  by  ~  20%  compared  to  the  same 
value  at  zero  temperature. 

Below  are  tabulated  results  of  computational  experiment  on  the  calculation 
of  the  theoretical  strength  for  TiB2  and  LaBe,  and  for  composite  based  on  them. 
Values  for  composite  are  calculated  using  the  rule  of  “mixtures”,  where  the 
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theoretical  strength  of  the  system  is  a  sum  of  the  theoretical  strengths  of 
components  taking  into  account  their  volume  fractions 

(^KZ  =  S  a(^a+Sb-(^T-  (62) 

Here  -  maximum  strength  of  MeB2,  -  strength  of  LaB6  at  maximum 
deformation  of  MeB2. 

In  Table  11  strength  (at  T  =  0)  corresponds  to  the  strength  of  LaBg  at 
deformation  of  the  crystal  up  to  about  0.1,  which  represents  the  ultimate  strain  for 
MeB2  (Me  -  Ti,  Zr,  Hf),  when  strength  of  diborides  reaches  maximum  values. 

As  a  result  of  the  calculation  using  the  formulas  (61,  62)  of  the  theoretical 
strength  of  the  components  and  composite  at  different  temperatures,  the 
dependence  on  the  temperature  of  the  theoretical  strength  under  uniaxial  tension 
for  LaBa,  TiB2,  as  well  as  LaB6-TiB2  system  was  obtained.  The  results  of  these 
calculations  are  presented  in  Table  1 1  and  in  Fig.  44  -  46. 

Table  11.  Temperature  dependence  of  the  theoretical  strength  (GPa)  of  LaB6 , 

TiB2,  as  well  as  LaB6-TiB2  composite. 


T,  K 

—  max 

0 

30.083 

47.180 

31.910 

300 

29.635 

46.804 

31.471 

500 

29.328 

46.639 

31.180 

750 

28.950 

46.271 

30.803 

1000 

28.610 

45.899 

30.460 

1500 

27.904 

45.230 

29.758 

2000 

27.461 

44.549 

29.290 

2500 

26.503 

43.672 

28.340 

2750 

26.203 

43.250 

28.027 

Theoretical  strength  of  LaBg  changes  almost  smoothly  with  the  temperature 
except  the  temperature  range  1300  K  <  T  <2200  K,  where  there  is  an  abnormal  part 
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of  dependence.  In  the  case  of  TiB2  such  interval  is  absent.  The  appearance  of 
anomalies  in  the  plot  of  the  theoretical  strength  from  the  temperature  is 
characteristic  for  LaBa,  and  it  affects  such  dependence  of  a  composite  LaB6-TiB2, 
\vhere  LaB6  content  is  75%  (mol.)  of  the  composite. 

It  should  be  noted  that  up  to  the  melting  temperature  borides  (TiB2,  ZrB2, 
HfB2)  and  composites  based  on  them  have  high  strength,  and  it  is  confirmed 
experimentally  by  S.  Ordanyan  [30]. 

If  the  crystal  is  subjected  to  deformation  with  simultaneous  raising  of  the 
temperature,  then  concurrently  with  deformation  a  mechanism  of  thermal 
expansion  will  also  operate.  Quantum  mechanical  calculations  show  that  the 
interaction  (binding)  between  the  boron  atoms  is  much  stronger  than  that  between 
metal  and  boron,  and  even  more  than  the  bond  of  metal  to  metal  [37  -  39],  ie,  to 
expand  B6  complex  it  requires  higher  temperatures.  With  thermal  expansion  and 
simultaneous  stretching  in  LaB6  there  is  a  change  of  the  relative  distance  between 
the  atoms  constituting  the  crystal,  which  in  turn  changes  the  characteristic 
temperature  range  where  there  is  an  abrupt  change  of  LCTE.  At  temperatures  T  > 
1300  K  the  expansion  of  Bg  complex  in  the  plane  [200]  resists  to  the  uniaxial  strain 
(stretching)  of  the  complex  along  the  perpendicular  direction,  with  the  result  that 
there  is  a  bump  on  the  plot  “theoretical  strength  -  temperature”. 

The  energy  of  the  electron-ion  system  (per  volume  of  the  unit  cell)  is  a 
function  of  the  crystal  temperature.  With  increasing  temperature,  energy  increases, 
and  the  maximum  strength  decreases.  We  introduce  the  concept  of  rate  of  strength 
change  with  temperature 

_3(7z 

.  (63) 

The  nature  of  the  dependence  of  j,  on  deformation  is  same  for  both  LaBg  and 

MeB2  (the  difference  is  only  in  parameters  -  deformations  and  the  corresponding 
strength)  (Fig.  47). 
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On  the  interval  (300  K  <  7  <  1300  K)  decreases  rapidly,  while 

transition  to  higher  temperatures  {T  >  1300  K)  change  in  deformation  does  not 
change  . 

(7  ,  GPa 


Fig.  44.  Dependence  of  theoretical  strength  of  LaB6  from  temperature  at  uniaxial 
deformation  (tension). 


GPa 


Fig.45.  Dependence  of  theoretical  strength  of  TiB2  from  temperature  at  uniaxial 
deformation  (tension). 
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Fig.  46.  Dependence  of  theoretical  strength  of  LaB6  -  TiB2  composite  from 
temperature  at  uniaxial  deformation  (tension). 
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Fig.  47.  Dependence  of  strength  change  rate  with  temperature  from  deformation  at 
uniaxial  tension  for  TiB2. 


Theoretical  strength  of  boride-boride  composites  and  their  constituent  parts 
varies  exponentially  as  a  function  of  temperature.  In  the  temperature  interval 
(1300  K  <  r  <2200  K)  for  LaB6  and  LaB6-TiB2  composite  the  observed  deviation 
in  the  dependence  “theoretical  strength  -  temperature”  is  a  consequence  of  changes 
in  the  coefficients  of  thermal  expansion  of  the  components  in  the  system. 
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Intercomponent  boundaries  and  energy  of  the  components  interface  in  quasi¬ 
binary  eutectic  systems 

A  study  of  interatomic  interactions  at  the  interface  of  eutectic  phases  was 
performed  in  [59,  60].  It  has  been  established  that  the  energy  states  of  atoms  in  the 
boundary  layers  and  volume  of  eutectic  phases  differ  considerably,  and  its  own 
interatomic  bonds  in  these  layers  are  also  changed.  These  data  confirm  that  the 
eutectic  does  not  represent  a  mechanical  mixture  of  phases,  but  a  system  of 
interacting  phases.  Regular  mutual  orientation  of  eutectic  phases  is  associated  with 
the  condition  of  minimum  interfacial  energy.  The  minimum  of  interfacial  energy  is 
achieved  when  the  combination  of  symmetry,  orientation  and  crystal  lattice 
parameters  of  the  two  crystals  corresponds  to  the  condition  of  the  largest  number 
of  atomic  coincidences  [61]. 

Under  the  surface  section  we  understand  a  transition  zone  which  separates 
from  each  other  two  phases,  or  (in  the  case  of  grain  boundaries),  two  crystals  of  the 
same  phase  [46]. 

Eutectic  solidified  under  normal  conditions,  has  a  surface  area  of 
approximately  1  m^  per  1  cm^  of  material  [46].  Extremely  large  interfacial  area, 
which  arose  in  the  process  of  crystallization,  is  a  characteristic  feature  of  eutectic 
compositions  and  determines  their  special  properties. 

To  determine  the  characteristics  of  the  interfaces  different  theories  are  used: 
1)  macroscopic,  such  as  thermodynamics,  allowing  us  to  obtain  the  averaged  data 
of  the  interface;  2)  atomic,  where  the  atoms  are  considered  as  spheres  with  a  short- 
range  (central)  forces  that  allow  to  draw  conclusions  about  the  structure  of  the 
interface  (to  build  geometric  and  energy  models);  3)  electronic,  which  give 
quantum-mechanical  description  of  the  interaction  between  the  atoms  at  the 
interface  in  the  eutectic  [46,  62]. 

A  comprehensive  description  of  all  of  the  observed  surface  phenomena  in 
eutectic  composites  at  the  interface  is  currently  not  available.  Therefore,  to  provide 
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a  general  overview  of  this  issue  is  difficult.  Furthermore,  most  studies  address 
grain  boundaries  or  external  surfaces  [46,  62];  phase  boundaries,  so  important  for 
the  eutectic  materials  are  just  beginning  to  attract  the  attention  of  researchers.  As 
the  authors  of  [62]  note,  the  experimental  data  for  heterogeneous  interface  are  not 
available.  For  this  reason,  the  description  of  the  interaction  between  the  atoms  of 
the  interface  is  a  complex  task,  especially  between  two  dissimilar  phases. 

The  interface  between  the  two  phases  can  be  non-coherent  or  coherent 
locally.  The  authors  of  [63,  64]  for  the  solution  of  this  problem  present  coherent 
local  areas  as  "cells".  For  a  complete  description  of  the  real  interface,  the  notion  of 
"supercell"  composed  of  the  crystal  lattices  of  the  two  phases  is  introduced. 

Based  on  the  results  of  calculations  from  first  principles  [37,  38]  for  metal- 
ceramic  composites  for  calculation  of  the  energy  of  the  electron-ion  system  of 
components  and  composite  as  a  whole,  a  model  was  created  to  describe  the 
formation  of  intercomponent  borders  in  the  quasi-binary  eutectic  systems. 

In  the  system  LaB6-MeB2  at  the  border  of  coupling  of  two  components  the 
binding  is  achieved  with  the  help  of  boron  atoms  belonging  to  two  components.  In 
the  case  of  LaB6-TiB2  composite  with  eutectic  composition,  there  is  0.75  part  of 
LaBe  molecule  and  0.25  part  of  TiB2  molecule  at  the  interface.  For  all  of  the 
system  as  a  whole  the  number  of  boron  atoms,  which  are  common  for  two 
components  at  the  border,  is  5.  Due  to  common  atoms,  the  boundary  is  strong 
enough.  The  total  number  of  these  common  atoms  depends  on  the  lattice 
parameters  of  the  components.  For  TiB2,  ZrB2  and  HfB2  these  parameters  are  very 
close,  and  so,  respectively,  energies  of  the  contact  are  nearly  identical,  which  is  not 
the  case  for  VB2  and  CrB2.  Here,  because  of  the  small  parameter  (c),  the  degree  of 
possible  coupling  of  two  components  increases  as  increases  the  probability  of 
coincidence  of  atomic  planes  [002],  [004]  LaBe  and  [002]  MeB2,  and  this  leads  to 
an  increase  in  energy  of  the  contact.  The  confirmation  of  these  statements  is 
obtained  from  the  results  of  calculations. 
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Based  on  the  fact  that  the  interface  is  associated  with  the  energy  of 
interaction  between  the  molecules  (atoms)  of  two  phases,  the  surface  energy 
associated  with  the  interface  of  two  phases  is  excess  internal  energy  of  the 
composite  compared  to  the  energy  of  the  components.  The  energy  of  the  phases 
(components)  and  the  system  based  on  them  are  calculated  using  the 
pseudopotential  method  [33,  41]. 

The  energy  of  components  calculated  per  one  molecule,  in  the  framework  of 
the  pseudopotential  method  is  obtained  by  adding  the  pair  intermolecular 
interactions,  which  are  determined  by  the  formula  [50] 


^{r)  =  U^^+^]^{q)Um{qr)qdq  , 
n  i  r 


where  Uec  -  electrostatic  energy,  and  O(^)  —  characteristic  function 

iq)  =  y^iq)%iq)£iq)  ■ 


(64) 


Here  V{q)  -  pseudopotential  of  the  system  (molecule)  [33],  e{q)  and  %{q)  - 
describe  screening  and  correlation  of  free  electrons  in  the  systems  under  study 
[50],  and  Q  -  volume  per  molecule.  Calculation  of  the  energy  component  was 
carried  out  under  the  scheme  [33],  namely,  the  total  energy  of  the  two  non¬ 
interacting  components  A  and  B  is  given  by 

U,=(C^-  U^  +  C-UJ+((l-Cf-U,,+C-U,)  ,  (65) 

where  include  kinetic  energy  of  free  electron  gas,  exchange-correlation 

effects,  energy  of  electrons  in  the  first-order  perturbation  theory  by  the 
pseudopotential  for  the  corresponding  molecules  [33];  -  energy  of 

interaction  between  A-A  and  B-B,  C  -  concentration  of  A  component. 

Energy  of  the  system  (A-B),  which  takes  into  account  interaction  between 
the  components  A  and  B  (Uab)  can  be  written  as 
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U,  =  C^-  U^  +  C-U,  +(l-Cf-U,,  +C-U,  +2-C-{\-C)-U^,. 


(66) 


Ultimately,  the  excess  energy,  which  determines  the  surface  factor,  is 

^U  ^U,-U,^-2-C\\-C)-U^,  .  (67) 

The  energy  of  interaction  between  the  two  components  {Uab)  can  be 
represented  as  the  sum  of  paired  intermolecular  potentials  (^,^  ) 

(68) 

where  Rij=  {A)-R.{B)  -  distance  between  the  molecules  i  -  j  belonging  to  the 

components  A  and  B.  Since  the  energy  of  interaction  between  two  molecules  does 
not  depend  on  the  position  of  other  molecules,  and  is  only  a  function  of  kp  (Fermi 
momentum)  or  Z/Q  [50],  it  is  possible  to  introduce  the  concept  of  "virtual"  cell 
with  a  volume 

Q=C-f2^  +  (l-C)-f25  (69) 

and  with  a  charge 

Z  =  C-Z,  +  (1-C)-Z^  ,  (70) 

where  Q,  {i  =  A  \B)  -  component  volumes,  Za\  Zb  -  respectively,  number  of  their 
outer  electrons.  Pseudopotential  of  the  virtual  crystal  will  be  equal  to 
U  =C  F^ +(1-C)  We  are  using  the  approximation  that  does  not  take  into 
account  variations  ;  Vb-'^v  the  nodes  of  the  unit  cell.  Under  this 

approach  the  characteristic  function  of  the  virtual  crystal  will  be 

^l{q)  =  (yv{q)f-lv{qHv{q)  •  (71) 

Functions  ey{q)  and  Xviq)  appearing  here  describe  screening  and  correlation  of 
electrons  for  the  virtual  crystal  and  are  determined  through  the  unit  cell  volume 
and  the  number  of  free  electrons. 
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Interaction  between  molecules  A  and  B  is  described  by  the  function 


KARv)  = 


(■Z,+-Z,y  ,  2Q, 


R. 


(271)- 


jo^^  (q)-&xp(iqRy)dq  ^  (72) 


where  Ry  -  intermolecula  distance  in  virtual  crystal. 

In  LaB6-MeB2  (Me  -  Ti,  Zr,  Hf)  systems  with  eutectic  content  the  content  of 
diborides  constitutes  a  small  percent.  We  assign  virtual  crystal  a  cubic  lattice.  Such 
unit  cell  contains  a  single  molecule,  which  consists  of  two  parts  -  molecule  of 
component  A  with  a  concentration  C  and  mplecule  of  component  B  with  a 
concentration  (1  -  C). 


Fig.  48.  In  the  unit  cell  of  the  virtual  crystal  the  surface  (MNPQ)  is  the  interface 
between  the  two  components 


Figure  48  presents  the  cube  (a  unit  cell  of  virtual  crystal)  and  the  surface  of 
coupling  of  two  components.  The  mean  value  of  the  contact  surface  area  in  a 
virtual  cell  can  be  determined  if  we  represent  Q  =  (virtual  cell  parameter  A  = 
MQ),  and  then  coupling  surface  area 

^  ~^MNPQ  ~  ^  •  (73) 

The  surface  energy  of  coupling  of  two  components  per  unit  area  of  contact  will  be 

y  =  -V^JS.  (74) 

Thus,  the  surface  energy  of  the  contact  between  two  insoluble  components 
will  be 

r=  -u„(R)ia^\  (75) 
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At  the  boundary  of  coupling  of  two  components  in  the  model  of  "virtual" 
crystal,  the  unit  cell  is  presented  as  averaged.  It  should  be  noted  that  in  the  quasi¬ 
binary  eutectic  systems  two  phases  or  two  components  (if  they  are  insoluble)  have 
different  crystal  lattice  as  well  as  sizes  of  atoms  or  molecules.  Accounting  for  these 
effects  in  the  calculation  of  the  surface  energy  of  interface  between  two 
components  in  the  framework  of  quantum-mechanical  calculations  is  quite 
challenging,  so  in  the  first  approximation  a  facet  of  virtual  cubic  unit  cell  of  the 
crystal  is  used  as  a  coupling  area  (Fig.  48). 

Table  12.  Calculated  volume  (Q),  contact  surface  (5”),  energy  of  interaction 
between  molecules  {Uab)  and  contact  energy  (y)  in  the  studied  eutectic  systems 


System 

Q,  a.u. 

S,  a.u. 

Uab,  a.u. 

y,  jW 

LaB6  -  TiB2 

406.107 

54.84 

0.12165 

3.454 

LaB6  -  ZrB2 

400.980 

54.12 

0.12420 

3.582 

410.786 

55.26 

0.11965 

3.388 

290.155 

43.83 

0.22886 

8.131 

LaB6  -  CrB2 

338.530 

48.57 

0.22516 

7.219 

B4C  —  TiB2 

576.077 

69.24 

0.80700 

18.149 

B4C-  SiC 

517.210 

64.50 

0.77800 

18.783 

TiB2-SiC 

151.326 

28.97 

1.58000 

40.050 

As  can  be  seen  from  the  obtained  results  (Table  12)  contact  surface  energy 
in  LaB6-MeB2  systems  is  small  compared  with  the  same  energy  in  B4C-TiB2,  B4C- 
SiC,  TiB2-SiC  systems.  For  LaB6-MeB2  (Me  -  Ti,  Zr,  Hf)  systems  contact  energy 
of  two  components  is  so  weak  that  can  be  ignored  and  the  composite  considered  as 
a  mechanical  mixture,  for  LaB6-VB2  and  LaB6-CrB2  systems  contact  energy  is 
nearly  2-2.5  times  higher.  In  other  cases,  such  as  SiC-B4C  binding  at  the  boundary 
is  realized  through  common  carbon  atoms.  B4C-  TiB2  system  can  have  shared 
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boron  atoms;  also  is  possible  formation  of  TiC  here,  but  this  is  unlikely  in  the 
presence  of  a  large  number  of  boron  atoms  belonging  to  two  components.  In  TiB2- 
SiC  system  the  bond  between  the  components  is  achieved  by  means  of  the 
emergence  of  TiC  clusters  locally  at  the  coupling  boundary,  and  in  this  case  the 
energy  of  the  contact  is  sufficiently  large,  and  this  is  confirmed  by  calculations. 

The  question  arises  —  is  it  characteristic  only  of  the  eutectic  composition,  or 
not?  One  can  definitely  answer  —  no.  For  these  systems  (at  least)  the  boundaries 
of  with  such  structure  will  take  place  regardless  of  the  composition.  Only  for 
eutectic  composition,  selected  from  the  condition  of  extremum  of  thermodynamic 
function,  contact  boundary  energy  minimizes  energy  of  the  system. 

Table  13.  Calculated  values  of  average  volume  and  of  contact  surface  energy  at 
different  temperatures  in  the  case  of  virtual  crystal  approximation 


LaB6-TiB2 

LaB6-ZrB2 

r,  K 

Q,  a.u. 

y,  J/m^ 

Q,  a.u. 

y,  J/m^ 

0 

406.107 

3.454 

400.980 

3.582 

300 

408.941 

3.407 

403.770 

3.535 

500 

410.830 

3.371 

405.613 

3.497 

750 

413.191 

3.331 

407.891 

3.457 

1000 

415.670 

3.290 

410.222 

3.416 

1500 

420.273 

3.214 

414.832 

3.336 

2000 

424.955 

3.139 

419.440 

3.271 

2500 

429.716 

3.056 

424.051 

3.109 

2750 

435.011 

3.010 

428.360 

3.106 

Using  the  model  proposed  in  [41]  (for  calculation  of  the  energy  of  electron- 
ion  system  of  the  crystal  in  the  second-order  of  perturbation  theory  with  respect  to 
pseudopotential  with  the  account  of  expansion  of  the  material  at  elevated 
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temperatures),  the  surface  energy  of  contact  was  calculated  depending  on  the 
temperature  for  the  temperature  range  from  zero  to  the  eutectic.  For  each 
temperature  the  energy  of  interaction  of  components  was  obtained  as  well  as  the 
volume  per  one  averaged  molecule,  and  with  the  help  of  (81),  (88)  the  energy  of 
the  contact  surface  at  different  temperatures  was  calculated.  Results  of 
computational  experiment  for  LaB6-TiB2  and  LaB6-ZrB2  systems  are  presented  in 
Table  13. 

With  increasing  temperature,  the  contact  surface  energy  decreases,  but  the 
system  retains  its  characteristic  mechanical  properties  up  to  the  melting 
temperature  (eutectic  temperature). 

So,  quasibinary  eutectic  is  a  system  consisting  of  two  almost  insoluble 
components,  bonded  on  at  the  boundary  with  the  help  of  a)  common  atoms  or  b) 
compounds  formed  from  atoms  of  the  components. 

The  energy  of  the  contact  surface  depends  on  ratio  between  the  parameters 
of  crystal  lattices  of  components. 

The  proposed  method  can  be  used  for  calculation  of  the  contact  surface 
energy  for  any  composite,  provided  components  are  insoluble. 


Accounting  for  the  effects  of  intercomponent  boundaries  in  LaB6-MeB2  (Me  -  Ti, 
Zr,  Hf)  systems  in  the  calculation  of  the  mechanical  properties 

The  idea  of  a  eutectic  as  a  simple  mechanical  mixture  of  fine  crystals,  gave 
reason  to  believe  that  mechanical  and  some  other  properties  of  alloys  related  by 
composition  to  the  two-phase  region  of  eutectic  systems  should  change  additively 
with  the  concentration.  According  to  the  well-known  scheme  proposed  by  N.S. 
Kumakov  and  S.F.  Zhemchzhnyi,  a  straight  line  should  represent  the  isotherm 
“structure-mechanical  properties”  of  these  alloys,  the  eutectic  point  for  such 
isotherm  does  not  stand  out  in  any  way.  However,  many  experimental  data 
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contradict  this  scheme.  The  common  rule  is  the  nonlinear  dependence  of  the 
mechanical  properties  from  the  composition  of  the  two-phase  alloy,  and  the 
additivity  may  occur  only  under  certain  conditions. 

In  [59,  60,  65]  a  number  of  studies  of  interatomic  interactions  at  the  interface 
of  eutectic  phases  was  carried  out.  It  was  established  that  the  energy  states  of  the 
atoms  in  the  boundary  layers  and  in  the  volume  of  eutectic  phases  differ 
considerably,  and  its  own  interatomic  bonds  in  these  layers  are  altered.  These  data 
confirm  that  the  eutectic  does  not  represent  a  mechanical  mixture  of  phases,  rather 
single  system  of  interacting  phases. 

Quasi-binary  eutectic  composite  is  a  system  consisting  of  two  pure 
components  (if  there  is  no  solubility  of  the  components)  and  contact  boundary 
between  them,  which  may  be  coherent,  semicoherent  and  amorphous  depending  on 
the  ratio  between  the  parameters  of  their  crystal  lattices. 

If  the  general  principle  of  calculating  the  theoretical  strength  under  uniaxial 
tension  through  the  energy  of  the  system  for  the  pure  components  is  known,  such  a 
task  is  quite  uncertain  and  complex  for  composites  and  it  is  due  to  the  influence  of 
the  boundary  of  coupling  between  components. 

Method  for  determination  of  the  total  energy  of  the  alloy  per  molecule 
within  the  pseudopotential  method  framework  is  based  on  the  summation  of  paired 
intermolecular  interactions  [33]. 

If  C  is  the  concentration  of  component  A,  then  for  the  energy  of  the  system 
from  two  components  one  can  write 

U  =  +U  BE  (\-Cf  +2C(\-C) U^, .  (76) 

Here 

(77) 
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(78) 


i*i 

C/,,-2C(l-C)  =  ^X‘J’«Wi).  (79) 

i*j 

where  Oaa,  ^bb,  ^ab  -  potentials  of  interaction  between  the  components  [35]. 

It  is  assumed  that  both  crystals  are  ideal,  the  first  component  has  only 
bindings  A  -  A,  the  second  component  bindings  B  -  B,  and  at  the  interface  —  only 
bindings  A  -  B.  Relation  (89)  can  be  written  as 

U  =  +  U,,  (1  - C)^  +  2 C (1  - C) U,,  =  +  {\-CfUl,,  (80) 


r^e 


Uaa  =  Uaa  +  Uab 


1-C 


V‘bb=+U,,+U,,— 


1-C 


Stresses  along  the  z-axis  are  defined  by  the  relation  [39] 

^  1  ac,,(j) 

S.d  de^  ’ 


(81) 


(82) 


(83) 


where  -  relative  deformation,  and  Uu  is  Uaa  or  Ubb,  d  -  lattice  parameter  in  the 
direction  of  the  axis  of  deformation.  Si  -  area  of  an  atomic  plane  in  a  crystal  lattice, 
perpendicular  to  the  axis  of  deformation. 

To  calculate  the  physico-mechanical  properties  of  the  composite  it  is 
necessary  to  find  the  dependence  of  the  energy  of  electron-ion  system  from  the 
lattice  parameter  of  the  system.  A  common  lattice  parameter  is  necessary  for  the 
quasi-binary  eutectic  composite,  which  is  absent,  because  each  component 
crystallizes  in  its  own  crystal  structure.  To  solve  this  problem,  we  introduce  the 
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concept  of  "virtual"  cell  along  the  coupling  boundary.  The  unit  cell  volume,  which 
corresponds  to  the  composition  of  such  composite  Q=  +(l-C)Qg  (C  - 
concentration  of  one  of  the  components  of  the  quasi-binary  composite),  then  the 
system  has  a  common  lattice  parameter  and,  accordingly,  the  calculation  of  the 
theoretical  strength  for  such  systems  is  carried  out  as  usual. 

If  we  present  virtual  crystal  as  a  mixture  of  two  components,  and  take  into 
account  that  each  cell  corresponds  to  a  single  molecule,  then  it  is  possible  to 
calculate  total  energy  of  the  system  similar  to  substitution  alloys. 

We  calculate  the  theoretical  strength  under  uniaxial  tension  by  standard 
procedure,  assuming  that  the  axis  of  deformation  of  the  virtual  crystal  is  parallel  to 
the  axis  of  deformation  for  pure  components.  After  calculating  the  strength  under 
uniaxial  deformation  (deformation  direction  is  fixed,  it  is  the  same  as  the  direction 
of  the  axis  of  the  fiber),  we  separate  the  contribution  from  member  Uab  (the  part 
that  is  responsible  for  the  interaction  A-B).  As  a  result,  we  have  the  following 
summands: 

^  _  1  ^  _  1  dUssids).  ^  _  1  dUAsidy)^ 

d^-S^  de^  ’  ^  dg-Sg  de^  ’  dy-Sy  de^ 

Here  Oa,  <Jb  -  theoretical  strength  of  the  components  taking  into  account  their 
crystal  lattice,  cIa,  dg,  dy  —  parameters  of  the  unit  cell  of  A  and  B  components  and 
virtual  cell  in  the  direction  of  the  axis  of  deformation,  Sa,  Sb,  Sy-  areas  of  atomic 
planes  in  the  respective  cells  perpendicular  to  the  axis  of  deformation. 

Effective  interaction  between  the  elements  of  A-A  taking  into  account  the 
influence  of  the  boundary  A-B  can  be  represented  as  the  first  term  in  (80),  and 
the  second  term,  respectively  —  as  an  effective  interaction  between  B  -  B. 

The  corresponding  theoretical  strengths  will  be: 


= 


1  ^^aa(^a)| 

d^-S^  de^ 


1-C  1  dU^,(dy) 

C  dy-Sy  de  ^ 


+ 


1-C 

C 


AB 


(85) 
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(86) 


1  SUM) 


C  1  du^(4,)  C 

H - =a„+ — o, 

1-Cd^S^  de^  l-C 


By  the  rule  of  mixtures,  excluding  the  influence  of  coupling  boundaries  of  two 
components,  the  strength  of  quasi-binary  systems  will  be 


(87) 


and  taking  into  account  the  influence  of  the  coupling  boundaries  we  have 


(88) 


Where  ^b  -  volume  fractions  of  the  components  of  the  eutectic 

composition  in  the  system. 

More  research  was  conducted  on  a  system  of  LaB6-ZrB2,  because  other 
components  —  TiB2  and  HfB2  —  by  the  structure  and  properties  are  very  close  to 
ZrB2.  Calculation  of  maximum  strength  and  deformation  was  conducted  with  and 
without  considering  the  influence  of  the  coupling  boundaries  between  components 
in  the  system.  Below  in  tabular  and  graphical  form  are  given  results  of 
computational  experiments. 

Fig.  49  and  50  show  functional  dependence  of  the  strength  on  the 
deformation  for  components  and  composite,  0/(5max)s  -  strength  of  the  matrix  at  the 
maximum  deformation  of  diboride,  o,(5*niax)s  -  same  strength  but  with  the  account 
of  the  influence  of  boundaries.  The  choice  of  maximum  deformation  was  made 
based  on  the  stress-strain  diagram  for  each  component. 
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or  GPa 


Fig.  49.  Dependence  of  theoretical  strength  from  deformation  for  1)  ZrB2, 
2)  LaB6-ZrB2,  and  3)  LaB6,  calculated  by  (83),  (87) 

Oj  GPa 


Fig.  50.  Dependence  of  theoretical  strength  from  deformation  for  1)  ZrB2,  2) 
LaB6-ZrB2  and  3)  LaB6,  calculated  by  (85),  (86),  (88) 

By  comparing  the  results  on  theoretical  strength  of  LaB6-ZrB2  composite 
obtained  from  computer  simulation  by  applying  the  rule  of  "mechanical  mixture" 
(Figure  49),  and  in  the  case  when  influence  of  the  contact  boundary  of  two 
components  (interface)  is  taken  into  account  (Figure  50),  we  can  conclude: 

1.  In  the  first  case  (the  rule  of  "mechanical  mixture")  characteristic  dependence 
strength-deformation  of  the  composite  and  the  matrix  (LaB6)  are  almost 
identical. 

2.  In  the  second  case,  the  shape  of  dependence  of  strength  on  the  composite 
deformation  differs  from  that  characteristic  both  for  matrix  (LaBg)  and  for 
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reinforcing  phase  (ZrB2).  Accounting  for  the  effect  of  the  interface  in 
calculations  leads  to  higher  values  of  mechanical  properties  of  the 
composite,  increasing  its  theoretical  strength. 

Table  14.  Maximum  deformation  and  theoretical  strength  (GPa)  for  components 
and  composite  without  additional  term  (Smax,  and  with  it  (e^max,  <3*max) 


Smax 

C^max 

(Smax)B 

S  max 

^  max 

ZrB2 

0.099 

39.060 

0.116 

43.170 

LaBe 

0.160 

31.589 

27.150 

LaB6-ZrB2 

0.099 

29.075 

0.116 

32.209 

Maximum  deformation  (Smax)  is  0.099  for  ZrB2  and  0.16  for  LaB6  without 
additional  term,  and  taking  into  account  the  influence  of  the  coupling  boundaries  of 
two  components  the  maximum  deformation  increases  to  0.1 161  in  the  case  of  ZrB2 
and  decreases  for  LaB6. 

If  we  assume  that  maximum  strength  of  the  matrix  in  its  pure  form  is  31.589 
GPa,  it  is  clearly  seen  that  strength  of  the  composite  with  the  influence  of  the 
contact  boundaries  is  greater  than  maximum  strength  of  the  pure  component  A. 

Using  formulas  (83)  (88)  strength  and  corresponding  deformation  are 
calculated  as  a  function  of  the  concentration  of  components  in  LaB6-ZrB2  system. 
Table.  15  shows  the  results  of  calculations  for  maximum  strength  for  hypoeutectic, 
eutectic  and  hypereutectic  composition. 

Table  15.  Dependence  of  maximum  deformation  and  strength  of  composite  (GPa) 
from  concentration  of  components  (Cb  =1  -  C  —  concentration  of  ZrB2) 


Cb 

0.2 

0.25 

0.28 

0.3 

0.31 

0.32 

0.35 

0.5 

Smax 

0.1181 

0.1171 

0.1165 

0.1160 

0.1155 

0.1150 

0.1140 

0.1110 

^max 

31.457 

31.560 

31.860 

32.027 

31.999 

31.999 

32.450 

33.703 

Distribution  A:  Approved  for  pubiic  reiease;  distribution  is  uniimited. 


With  increasing  of  ZrB2  concentration  the  maximum  deformation  of  the 
composite  decreases  and  strength  increases.  When  Cb  =  0.3  the  composite  has 
maximum  strength  in  the  region  0  <  Cb  <  0.33.  A  further  increase  in  concentration 
leads  to  ZrB2  gradual  increase  in  the  strength  of  the  composite,  and  maximum 
value  is  reached  at  Cb  =  1. 

Account  of  the  interaction  of  components  in  the  composites  at  the  coupling 
boundary  leads  to  an  increase  of  both  the  maximum  deformation  of  refractory 
component  in  the  composite  and  the  theoretical  strength  of  the  composite. 

The  obtained  local  maximum  of  strength  for  the  composite  in  the  eutectic 
point  is  the  confirmation  of  the  fact  that  eutectic  is  not  just  a  mechanical  mixture 
and  that  for  the  eutectic  are  characteristic  not  only  minimum  temperature  of 
formation,  but  also  superiority  of  mechanical  properties. 

Conclusion: 

-  A  model  is  developed  to  account  for  changes  in  the  volume  of  the  unit  cell  at 
temperatures  different  from  zero  based  on  the  method  of  pseudopotentials 
(second-order  of  perturbation  theory)  the  so-called  quasi-harmonic 
approximation. 

-  -  In  quasi-harmonic  approximation  the  calculated  thermal  expansion 
coefficients  for  LaB6  and  LaB6-MeB2  (Me  -  Zr,  Ti)  systems  based  on  it  have  a 
jump  in  the  temperature  range  from  750  K  to  1000  K,  which  is  associated  with 
the  presence  of  a  deeper  minimum  of  energy  of  Be  complex  than  interaction 
energy  between  La  and  B  atoms,  or  La-La  atoms.  Therefore,  in  the  initial  stage 
of  raising  the  temperature  to  750  K  it  is  La  atoms  that  are  displaced  from  the 
equilibrium  position.  At  higher  temperatures  Be  volume  complex  begins  to 
participate  in  the  expansion,  which  leads  to  a  jump  in  the  coefficient  of  thermal 
expansion. 
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With  the  help  of  quasi-harmonic  approximation  theoretical  strength  of  borides 
(LaB6,  MeB2)  and  LaB6-TiB2  system  with  eutectic  composition  under  uniaxial 
deformation  in  the  temperature  range  from  0  to  2750  K  were  calculated.  The 
characteristic,  almost  linear  dependence  of  the  strength  from  the  temperature 
changes  in  the  temperature  interval  (1300  K-2200  K)  in  the  case  of  LaBe,  an 
increase  in  temperature  leads  to  an  increase  in  the  theoretical  strength.  Same 
dependence  is  observed  for  LaB6-MeB2  (Me  -  Ti,  Zr,  Hf)  systems,  where  in 
eutectic  systems  LaBe  is  0.7  mol.  part.  The  observed  deviation  in  the 
dependence  “theoretical  strength  -  temperature”  is  a  consequence  of  change  in 
the  coefficient  of  thermal  expansion. 

Based  on  the  results  of  computational  experiments,  a  model  of  formation  of 
interphase  boundaries  in  the  quasi-binary  eutectic  systems  was  proposed. 
Quasibinary  eutectic  is  a  system  consisting  of  two  almost  insoluble 
components,  bonded  on  the  coupling  boundary  using  shared  atoms  or 
compounds  formed  from  atoms  of  the  components. 

With  the  help  of  introduced  concept  of  "virtual  cell"  along  the  coupling 
boundary  of  components  the  surface  energy  of  the  interface  between  the 
components  in  the  boride  and  metal-ceramic  quasi-binary  composites  with 
eutectic  composition  was  calculated.  The  surface  energy  of  the  interface  for 
LaB6-TiB2,  LaB6-ZrB2  systems  at  different  temperatures  was  calculated.  With 
increasing  temperature,  the  contact  surface  energy  decreases,  but  the  system 
retains  its  characteristic  mechanical  properties  up  to  melting  temperature 
(eutectic  temperature). 

The  contact  surface  energy  depends  on  the  ratios  between  the  crystal  lattice 
parameters  of  components. 

A  relationship  was  obtained  that  takes  into  account  effect  of  the  interface  on 
mechanical  properties  of  the  composites.  Using  these  relations,  theoretical 
strength  of  LaB6-ZrB2  and  maximum  possible  deformation  were  calculated. 
Under  this  approach,  the  dependence  of  the  theoretical  strength  on 
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concentration  of  components  was  found,  and  the  existence  of  a  local  maximum 
of  theoretical  strength  for  the  composite  in  the  eutectic  point  was  demonstrated. 

-  The  obtained  local  maximum  of  strength  for  the  composite  in  the  eutectic  point 
is  the  confirmation  that  eutectic  is  not  just  a  mechanical  mixture  and  it  is 
characterized  by  not  only  the  minimal  temperature  of  formation,  but  also  the 
superiority  of  mechanical  properties. 


6.  The  ideal  work  of  adhesion  interface  component  eutectic 

One  of  the  promising  methods  used  in  materials  science  for  production  of 
composite  materials  is  the  directional  solidification  of  eutectic  ceramic. 
Directionally  crystallized  under  certain  process  conditions  composites  based  on  of 
LaB6-TixZri_xB2  system  have  a  high  level  of  physical  and  mechanical  properties 
such  as  strength  and  fracture  toughness.  This  level  of  properties  is  achieved  due  to 
specific  eutectic  microstructure  of  the  composite  —  diboride  fibers  uniformly 
distributed  in  the  matrix  of  lanthanum  hexaboride.  Computer  simulation  of  the 
formation  of  these  eutectic  alloys  using  phenomenological  models  allows  us  find 
optimal  modes  of  crystallization.  However,  these  models  include  a  number  of 
parameters  which  are  very  difficult,  and  sometimes  impossible,  to  establish 
experimentally.  One  such  parameter  is  the  ideal  work  of  adhesion,  which  is  defined 
as  the  work  required  for  reversible  separation  of  the  interface  into  two  free  surfaces 
without  plastic  deformation  and  diffusion.  From  a  practical  point  of  view,  of 
considerable  interest  is  to  clarify  the  nature  of  the  influence  of  impurities,  in 
particular  titanium  on  the  properties  of  the  interface.  It  was  established 
experimentally  that  between  the  components  of  LaB6-TixZri_xB2  system  the 
following  crystallographic  orientation  relationships  are  observed:  [0011- 
LaB6//[001]-TixZri_xB2,  (110)-LaB6//(110)-TixZri_xB2.  For  these  orientation 

relationships  the  results  of  calculations  are  presented  of  ideal  work  of  adhesion  for 
both  coherent  and  semicoherent  interfaces  for  v  =  0,  0.57,  1. 
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The  ideal  work  of  adhesion  Wad  is  defined  as  the  energy  required  to  separate 
interface  into  two  free  surfaces.  At  the  same  time  plastic  and  diffusion  processes 
are  neglected  [66].  In  accordance  with  this  definition,  adhesion  energy  can  be 
calculated  as  follows: 


W^,={E,  +  E,-E,,)IS,  (82) 

where  -  total  energy  of  supercell  containing  interface  LaB6  and  TixZri.xB2.  E^ 
and  Eg  -  total  energies  of  isolated  phases  LaBe  and  TixZri.xB2  respectively,  S  - 
contact  area  KOHxaKxa. 

The  calculations  of  the  total  energies  was  carried  out  using  the  software 
package  Quantun-Espresso  (version  4.3.2)  [67].  This  package  implements  ab  initio 
pseudopotential  method  in  combination  with  density  functional  theory.  To  obtain 
adequate  results  with  ab  initio  methods  we  need  carefully  select  a  number  of 
calculation  parameters  that  determine  the  accuracy  of  calculations.  To  such 
parameters  in  the  ab  initio  pseudopotential  method  belong  cutoff  energy  Ecut  in  the 
expansion  of  plane  waves  and  the  number  k  of  points  in  the  Brillouin  zone.  To 
select  these  parameters  calculations  of  the  total  energy  of  LaB6,  ZrB2  were  carried 
out.  Calculations  were  carried  out  for  experimental  parameter  data.  The  results  of 
calculations  are  presented  in  Fig.  51,  and  variants  of  systems  of  points  k  —  in 
Table  16. 

Analysis  of  the  results  shows  that  the  accuracy  of  the  calculation  of  the  total 
energy  1  mRy  is  achieved  by  using  a  system  of  points  16x16x16  (variant  5),  Ecut  = 
40  Ry  for  LaB6  and  Ry  =  30  Ry  for  ZrB2.  These  values  were  used  to  calculate  the 
equilibrium  lattice  parameters  of  LaBe,  ZrB2  and  TiB2. 
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Fig.  51.  Total  energy  of  LaB6  and  ZrB2  at  different  temperatures  of  calculations. 


Table  16.  Calculation  parameters. 


Calculation  variant 

System  of  points  k 

1 

6x6x6 

2 

8x8x8 

3 

10x10x10 

4 

12x12x12 

5 

16x16x16 

6 

20x20x20 

LaB6  belongs  to  a  simple  cubic  structure,  however,  it  has  two  structural 
parameters,  lattice  parameter  a  and  parameter  u,  which  determines  the  position  of 
boron  atoms  in  the  unit  cell.  For  their  determination  calculations  of  the  total  energy 
for  a  set  of  lattice  parameters  with  relaxation  of  the  positions  of  boron  atoms  were 
carried  out.  The  obtained  dependence  is  shown  in  Figure  52. 
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Fig.  52.  Dependence  of  total  energy  of  LaB6  from  lattice  parameter. 

The  equilibrium  lattice  parameter  determined  from  this  dependence  is  <2o(LaB6)  = 

o 

4.1491  A.  This  lattice  parameter  corresponds  to  the  value  u  =  0.1999.  The 

o 

experimental  values  are,  respectively,  0.1995  and  4.1491  A  [68]. 

ZrB2  and  TiB2  belong  to  hexagonal  structure  and  have  two  lattice  parameters 
a  and  c.  Total  energy  of  ZrB2  was  calculated  on  the  mesh  from  five  values  of  a  and 
five  values  of  da.  The  figure  shows  surface  contours  of  the  total  energy  in  the 
coordinates  volume  (V)  -  ratio  da.  The  resulting  surface  was  approximated  by 
two-dimensional  Chebyshev  polynomials  of  order  4.  The  accuracy  of  the 
approximation  is  2-4  units  of  sixth  decimal  place  in  the  total  energy.  Determination 
of  the  minimum  of  approximating  function  yielded  the  following  values  of  the 

o  o 

lattice  parameters:  ao^LxE^  =  3.1697  A,  co(ZrB2)  =  3.5418  A.  The  experimental 
values  are  respectively  3.170  and  3.532  A  [4].  For  TiB2  calculated  values  of  the 
lattice  parameters  were  <2o(TiB2)  =  3.0307  A,  co(TiB2)  =  3.2208  A  (experimental 
values  are  equal  to  a  =  3.026  A,  c  =  3.213A  [69]). 
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Fig.  53.  Contours  of  the  surface  of  total  energy  for  ZrB2. 


For  eutectic  of  LaB6-TixZri_xB2  the  following  orientation  relationships  in  the 
formation  of  the  interface  between  the  diboride  fibers  and  lanthanum  hexaboride 
matrix  were  experimentally  observed:  [001]-LaB6  //  [0011-ZrB2,  (110)-LaB6  // 
(110)-ZrB2.  A  detailed  analysis  of  these  orientation  relationships  gives  the 
following  results.  For  lanthanum  hexaboride  at  the  interface  plane  exists  a 
rectangular  lattice  of  lanthanum  atoms  with  parameters  <2(LaB6)  =  V2<2o(LaB6)  = 
5.8677,  Z7(LaB6)  =  <2o(LaB6)  =  4.1491  A.  Each  node  of  the  plane  lattice  is 
associated  with  boron  octahedron  with  two  boron  atoms  lying  in  the  plane  and  the 
other  four  are  located  symmetrically  on  opposite  sides  of  the  plane. 

[0011 


Fig.  54.  Interface  plane  of  LaB6. 
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For  zirconium  diboride/titanium  on  interface  plane  also  exists  a  rectangular  lattice 
of  atoms  zirconium/titanium  with  parameters  fl(ZrB2)  =  V3flo(ZrB2)  =  5.4901, 
b{ZxB2)  =  co(ZrB2)  =  3.5418  A  (fl(TiB2)  =  V3ao(TiB2)  =  5.2493,  Z7(TiB2)  =  co(TiB2) 
=  3.2208  A). 


[0011 


Thus  the  lattice  mismatch  at  the  interface  planes  of  LaB6  and  ZrB2/TiB2  turns  out 
to  be  different  in  the  directions  [TlO]  and  [0011.  So,  for  the  [TlO]  direction  this 
discrepancy  is  5i(ZrB2)  =  0.066  (5i(TiB2)  =  0.111),  and  for  the  [0011  direction 
52(ZrB2)  =  0.158  (52(TiB2)  =  0.252).  These  values  show  that  it  is  possible  for  ZrB2 
in  the  direction  of  [TlOl  provide  a  coherent  interface  with  a  relatively  small 
mismatch.  The  lowest  commensurate  unit  cell  along  the  [001]  direction  should 
consist  of  six  unit  cells  of  LaB6  and  seven  cells  of  ZrB2  and  contain  misfit 
dislocations.  For  TiB2  mismatch  in  both  directions  is  significantly  greater  than  for 
ZrB2.  Nevertheless,  both  coherent  and  semicoherent  interface  for  zirconium  and 
titanium  diboride  were  constructed  the  same  way. 

To  perform  calculations  of  ideal  work  of  adhesion  Wad  as  a  first  approximation 
we  developed  four  models  of  interface  M1-M4,  coherent  along  the  [TlO]  and  [001] 
directions  (Fig.  56).  In  order  to  ensure  coherency  the  lattice  parameters  of  ZrB2  and 
TiB2  were  changed  accordingly.  Supercell  consisted  of  three  layers  of  LaB6  and 

o 

three  layers  of  ZrB2/TiB2,  separated  by  a  vacuum  gap  of  ~  10  A  and  contained  70 
atoms. 
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O-Zr 


o-Ti 
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Fig.  56.  Atomistic  models  of  interface  [001]-LaB6//[001]-Zr(Ti)B2,  (110)- 
LaBe/ZCl  10)-Zr(Ti)B2:  a  -  d  —  models  Ml  -  M4  respectively 

In  the  Ml  model  boron  atoms  on  the  interface  surfaces  of  LaB6  and  ZrB2  are  linked 
by  “bridging”  ties.  Due  to  this  structure  of  interface,  there  are  regions  in  which 
atoms  of  lanthanum  and  zirconium  are  situated  opposite  each  other  and  are  located 
close  enough.  In  the  M2  model  lanthanum  and  zirconium  atoms  are  surrounded  by 
atoms  of  boron  and  do  not  contact  each  other.  In  the  M3  and  M4  models  zirconium 
atoms  closest  to  the  interface  plane  are  replaced  by  titanium  atoms.  This 
substitution  allows  us  to  trace  the  influence  of  titanium  admixture  on  the  ideal 
work  of  adhesion. 

The  energies  of  the  systems  under  consideration  are  calculated  as  a  function  of  the 
distance  between  the  interface  surfaces.  The  obtained  dependences  are  shown  in 
Fig.  57.  From  these  curves  the  equilibrium  distances  do  between  the  interface 
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surfaces  and  the  magnitude  of  ideal  work  of  adhesion  Wad  are  determined.  The 
results  of  the  calculation  are  given  in  Table  17. 


Fig.  57.  Interface  energy  as  a  function  of  the  distance  between  the  interface 
surfaces:  a  -  d  —  models  Ml  -  M4,  respectively. 

Table  17.  Ideal  work  of  adhesion  and  equilibrium  distances  between  interface 
planes. 


Parameter 

Ml 

M2 

M3 

M4 

1.89 

1.34 

2.23 

5.39 

cIq,  nm 

0.189 

0.198 

0.186 

0.127 

At  the  equilibrium  distance  do  =  0.189  nm  between  the  interface  surfaces,  the 
distance  between  the  nearest  atoms  of  zirconium  and  lanthanum  in  the  model  Ml  is 
equal  to  0.277  nm.  This  distance  is  less  than  the  sum  of  the  atomic  radii  (0.160  nm 
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for  Zr  and  0.186  nm  La  [70])  of  the  respective  elements.  In  connection  with  this, 
the  question  arises  of  the  validity  of  the  model  Ml.  Note,  however,  that  in  LaBe 
and  ZrB2  chemical  bond  between  the  metal  atoms  and  boron  has  a  partially  ionic 
character  and  the  metal  atoms  are  ionized  in  these  compounds.  It  is  therefore  more 
correct  to  consider  the  ionic  radii  which  are  equal  to  0.087  nm  for  Zr”^"^  and  0.122 
nm  for  La"^^  [70].  Comparing  Fig.  57a,  b  and  57c,  d  shows  that  more  energetically 
advantageous  are  interface  models  M4  and  M2  in  comparison  with  the  models  Ml 
and  M3,  respectively. 

From  these  data  it  follows  that  the  presence  of  titanium  atoms  in  the  plane  of 
the  interface  leads  to  a  considerable  increase  in  the  work  of  adhesion. 

The  ideal  work  of  adhesion  Wad  in  the  approximation  of  semicoherent 
interface  taking  into  account  relaxation  of  atomic  positions  were  calculated  for 
ZrB2,  Tio.57Zro.43B2,  TiB2. 

Fig.  58  shows  the  model  of  semicoherent  interface. 


O-La  o-Zr  •-B 


Fig.  58.  Model  of  semicoherent  interface  for  LaB6-ZrB2:  left  -  top  view;  right  - 
side  view. 
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For  calculations  of  ideal  work  of  adhesion  Wad  in  the  approximation  of 
semicoherent  interface  taking  into  account  relaxation  of  atomic  positions  in  a 
reasonable  time,  we  had  to  use  a  supercell  with  fewer  layers  of  LaB6,  ZrB2/TiB2, 
than  for  calculations  in  the  approximation  of  a  coherent  interface.  However,  the 
computational  experiments  showed  that  a  reduction  in  the  number  of  atoms  in 
supercell  by  30%  changes  only  third  decimal  place  in  the  calculated  ideal  work  of 
adhesion.  In  this  case,  the  calculation  time  is  reduced  by  about  ten  times 

Fig.  59  shows  the  appearance  of  semicoherent  interface  LaB6-ZrB2  before 
relaxation  and  after  relaxation  of  atomic  positions. 


o  o  o  o 


o  o  o  o 


o  o  o  o 

0  0  0  ° 


o  o  o  o 


o  O  o  O 


a)  b) 

Fig.  59.  Semicoherent  interface:  a)  before  relaxation  of  atomic  positions;  b)  after 
relaxation.  Atoms  La  -  large  gray  circles,  Zr  -  small  gray  circles,  B  -  black  circles. 

Energy  of  LaB6-TixZri_xB2  (x  =  0,  0.43,  1)  systems  was  calculated  as  a  function  of 
the  distance  d  between  the  interface  surfaces  (Fig.  60). 
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Fig.  60.  Energy  of  the  interface  as  a  function  of  the  distance  d  between  interface 
surfaces. 


The  table  below  shows  the  equilibrium  distances  do  between  interface  surfaces  and 
ideal  work  of  adhesion  Wad  calculated  from  these  curves. 


Table  18. 


X 

0 

dg  5  A 

0.00 

1.429 

3.828 

0.57 

1.329 

3.874 

1.00 

1.139 

4.075 

The  table  shows  that  the  presence  of  titanium  increases  the  ideal  job  of  adhesion. 
The  greatest  work  of  adhesion  is  observed  for  the  interface  with  TiB2. 

In  [71]  it  was  shown  that  the  substitution  of  zirconium  by  titanium  causes  a 
change  in  the  growth  direction  of  the  fibers  during  directional  solidification  of 
LaB6-TixZri_xB2  eutectics.  One  reason  for  this  behavior  of  the  system  may  be  a 
significant  energy  loss  in  case  of  constant  growth  direction  of  the  TiB2  fibers.  For  a 
qualitative  understanding  of  the  problem  the  calculations  have  been  carried  out  of 
the  magnitude  of  additional  energy  that  have  to  be  expended  to  ensure 
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semicoherent  interface  for  TiB2  and  ZrB2  taking  into  account  the  availability  of 
free  surface.  For  ZrB2  additional  energy  was  1.11  eV/(elem.  cell),  and  for  the  TiB2 
-3.78  eV/(elem.  cell).  Thus,  energy  loss  during  the  formation  of  the  semicoherent 
interface  in  the  case  of  titanium  diboride  is  substantially  higher  -  3.4  times  than  in 
the  case  of  zirconium  diboride. 

Conclusions. 

The  ideal  work  of  adhesion  of  interface  [001]-LaB6//[001]-TixZri_xB2,  (110)- 
LaB6//(110)-TixZri_xB2  increases  with  increasing  titanium  content.  Change  in  the 
direction  of  fibers  growth  during  directional  solidification  of  LaB6-TixZri_xB2 
eutectics  is  associated  with  a  large  loss  of  energy  in  the  formation  of  semicoherent 
interface  with  increasing  titanium  content. 

7.  Simulation  of  the  segmentation  process  of  Ti^CB  and  Tii3Ci4  nanoclusters 

The  interest  to  cluster  models  of  crystal  segmentation  increased  after  stable 
nanoclusters  were  experimentally  obtained.  For  the  Ti-C  system  TisC^ 
titanocarbohedrene  [72]  proved  to  be  the  most  stable,  having  the  fullerene-like 
form  of  C20.  It  was  also  shown  [73],  that  clusters  with  such  stoichiometric 
composition  can  exist  in  different  symmetric  forms.  In  [74]  the  synthesis  of  a 
number  of  cubic  nanocrystals  with  the  structure  of  NaCl  was  reported:  TiuCo 
(3x3x3),  TiigCis  (4x3x3),  Ti23C22  (5x3x3),  Ti24C24  (4x4x3),  Ti27C27  (6x3x3), 
TiaoCso  (5x4x3),  Ti32C32  (4x4x4),  Ti36C36  (6x4x3).  The  existence  of  these  samples 
was  confirmed  in  experiments  [75]. 

In  [76]  “ab  initio”  calculations  were  carried  out  within  the  framework  of 
density  functional  theory  using  the  method  of  generalized  gradient  approximation 
(GGA)  for  a  series  of  clusters  including  mentioned  Ti^Co,  Ti^Cu  and  truncated 
fragments  of  these  clusters  without  comer  atoms  Ti6Ci3  and  Ti^Ce.  With  the  help 
of  estimates  of  the  interatomic  distances  for  Ti-C  bonds  and  average  energies  of 
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these  bonds  it  was  concluded  that  angle  atoms  weakened  by  bonds  do  TiiaCu 
cluster  metastable.  This  conclusion  is  also  confirmed  by  calculations  of  the 
vibrational  spectra  of  atoms. 

Obviously,  the  combination  of  many  unique  structural  and  electronic 
features,  it  creates  the  preconditions  for  the  development  of  materials  promising 
for  application  in  many  fields  by  segmenting  or  assembling  clusters  [771.  Self- 
assembly  of  nanoclusters  in  crystalline  structures  strongly  depend  on  the  type  of 
symmetry  of  nanoclusters  themselves  and  composites  based  on  them.  In  [781 
group-theoretic  relations  are  detailed  conducive  to  the  evolution  of  the  crystals.  It 
is  shown  that  the  same  group-theoretic  conditions  can  lead  to  a  variety  of 
macrostructures,  determined  by  the  physico-chemical  properties  of  atoms  (e.g. 
atomic  radii,  electronic  structure),  and  this  stipulated  the  introduction  of  the 
concept  of  channels  of  evolution. 

Questions  of  self-assembly  or  segmentation  of  nanocrystals  mentioned 
above,  as  well  as  their  stability  has  not  yet  been  studied  adequately.  In  particular,  it 
is  interesting  how  cubic  nanocrystallites  of  (3x3x3)  (4x3x3)  type  and  the  like,  but 
with  a  larger  number  of  atoms,  which  are  composites  of  a  series  of  two  flat 
fragments  of  Ti5C4  and  Ti4C5  are  formed.  Can  there  be  such  forms,  at  least  in  a 
metastable  state;  whether  fragmentation  of  nanocrystals  from  these  planes  is 
possible?  In  this  section  we  describe  approaches  to  answer  these  questions.  With 
this  purpose  in  mind,  we  theoretically  investigated  the  stability  and  the  interaction 
of  plane  nanoclusters  of  Ti5C4  and  T^Cs  types  and  conducted  a  number  of 
computational  experiments  to  help  understand  the  conditions  of  their  segmentation 
into  cubic  clusters  TiuCia  and  Tii^Cu.  Cluster  model  calculations  were  carried  out 
“from  first  principles”  within  the  density  functional  theory  framework  using  the 
software  package  “Gaussian'03”  [791. 
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Electronic  properties  and  stability  of  planar  nanoclusters  Ti5C4  and  T^Cs 

Fig.  61  presents  planar  clusters,  which  can  be  regarded  as  the  splitting  of  the 
cubic  nanocrystallite  (3x3x3)  into  three  planes  Ti5C4-Ti4C5-Ti5C4  or  Ti4C5-Ti5C4- 
T^Cs. 

9  8  7 

6  5  4 

3  2  1 

a  -  Ti5C4;  b  -  T^Cs 

Fig.  61.  Optimized  planar  clusters 

In  the  first  case  TiuCn  crystallite  contains  one  atom  of  titanium  in  excess, 
and  in  the  second  Ti^Cu  -  one  carbon  atom.  As  can  be  seen  from  Fig.  61,  both 
planar  nanoclusters  have  only  three  inequivalent  atoms:  one  central  with  a 
coordination  number  K  =  4,  four  comer  (K  =  2)  and  four  atoms  at  the  centers  of  the 
edges  (K  =  3).  Accordingly,  they  form  12  bonds  Ti-C,  and  among  them  only  two 
are  different  “corner  atom  -  center  of  the  edge”  (1,2)  and  the  “center  of  the  edge  - 
central  atom  of  the  facet”  (2,5).  The  angle  between  the  atoms  2  -  1  -  4  we  denote 
as  A(2,l,4).  Table  19  provides  the  stmctural,  energetic,  and  electronic  properties  of 
these  forms.  The  average  energy  of  the  Ti-C  bond  was  calculated  as  the  binding 
energy  divided  by  the  number  of  bonds  in  the  cluster.  The  gap  in  the  energy 
spectmm,  which  separates  the  occupied  states  from  free  (LUMO-HOMO)  has  a 
plus  sign  in  the  presence  of  a  semiconducting  gap  and  the  minus  sign  in  the  case  of 
overlapping  zones.  The  values  of  the  electronic  charge  transfer  according  to 
Mulliken  are  given  for  non-equivalent  atoms. 
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Table  19.  Geometric  and  electronic  properties  of  planar  nanoclusters 


Ti5C4 

T^Cs 

Ti-C  Bond  length  (1.2),  nm 

0.1936 

0.1925 

Ti-C  Bond  length  Ti-C  (5.2),  nm 

0.2018 

0.1977 

Angle  (2,1,4)  grad. 

94.99 

93.13 

Binding  energy/atom,  eV 

8.3753 

8.0917 

Average  energy  of  Ti-C  bond,  eV 

6.2815 

6.0688 

LUMO-HOMO,  eV 

0.6942 

0.6095 

Charge  transfer  by  Mulliken: 

Angle  atom  (1) 

Ti  0.2117 

C  -0.2284 

Atom  in  the  edge  center  (2) 

C  -0.3135 

Ti  0.4071 

Ti  0.2589 
C  -0.3103 

Central  atom  (5) 

The  effect  of  optimization  can  be  seen  from  the  geometrical  parameters  of 
the  planar  clusters.  As  the  initial  value  of  the  Ti-C  bond  length  the  magnitude  of 
0.21  nm  was  taken  (the  experimental  value  0.2176  nm).  Table  19  shows  that  in  the 
result  of  optimization  the  compression  of  the  structure  has  occurred,  but  not 
entirely  uniform.  The  distance  from  the  central  atom  to  the  midpoints  of  the  sides 
of  the  rectangle  became  greater  than  the  distance  from  the  corner  atoms  towards 
the  midpoints  of  the  edges.  For  Ti5C4  the  ratio  of  these  lengths  is  equal  to  1.04,  and 
for  Ti4C5  -  slightly  less  than  1.03.  Accordingly,  the  angle  of  the  rectangle  in  the 
first  case  increased  to  95°,  and  in  the  second  -  to  93°.  Overall,  in  both  cases  the 
structure  of  the  figure  is  preserved  but  only  with  some  distortion.  Cohesive  energy 
and  the  average  energy  of  the  Ti-C  bond  is  a  little  more  in  the  cluster  Ti5C4.  This 
indicates  that  the  cluster  has  a  more  strong  structure.  Calculations  show  that  both 
objects  have  in  the  energy  spectrum  a  small  semiconductor  gap.  Charge  transfer  is 
of  the  same  type:  titanium  atoms  give  some  fraction  of  the  electron  charge  to 
carbon  atoms.  However,  if  in  Ti5C4  titanium  atom  located  in  the  center  of  the 
fragment,  gives  the  greatest  charge,  and  angular  titanium  atoms  -  lower,  but  equal. 
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and  carbon  atoms  have  equal  charge,  then  in  Ti4C5  titanium  atoms  4  and  6  give 
somewhat  greater  charge  than  atoms  2  and  8.  In  the  carbon  subsystem  central  atom 
takes  the  largest  share  of  the  charge,  and  angular  atoms  -  a  smaller  but  equal  parts. 

Software  package  “Gaussian'03”  allows  to  calculate  vibrational  spectra  of 
atoms,  molecules  and  clusters.  Calculation  of  atomic  vibrations  allows  one  not 
only  determine  the  most  active  frequencies  and  type  of  oscillation,  but  also  provide 
information  on  the  structural  stability  of  the  object.  Thus,  the  presence  of 
imaginary  frequencies  in  the  spectrum  indicates  that  during  optimization  not 
stationary  condition  was  found  but  transit  [801.  In  the  case  of  planar  nanoclusters 
the  calculations  indicated  the  presence  of  imaginary  frequencies  in  the  spectra  of 
both  structures.  It  means  that  Ti5C4  and  Ti4C5  are  metastable  structural  formations 
that  have  a  limited  lifetime.  During  this  period  of  time  nanoclusters  may  undergo 
structural  changes  or  merge  with  other  clusters  during  segmentation.  Oscillations 
corresponding  to  imaginary  frequencies  in  both  cases  are  due  to  oscillations  of 
carbon  subsystem.  In  the  case  of  Ti5C4  it  is  double  type  Eu  of  longitudinal  and 
transverse  (Z  =  0)  plane  oscillation  of  carbon  atoms.  In  Ti4C5  imaginary  frequency 
corresponds  to  oscillations  A2U  of  five  carbon  atoms  along  the  axis  Z, 
perpendicular  to  the  plane  of  the  nanocluster. 

Segmentation  of  planar  nanoclusters  Ti5C4  and  TUCs 

Geometric  optimization  of  a  molecule  or  cluster  in  software  package 
“Gaussian'03”  is  achieved  by  finding  the  minimum  of  the  potential  energy  surface 
by  gradient  methods  [801.  As  the  minimization  parameters  we  used  nonequivalent 
nearest  interatomic  distances,  angles  between  the  interatomic  bonds  and  sums  of 
the  lengths  of  the  bonds  between  the  four  atoms.  Thus,  even  for  small  clusters  the 
number  of  parameters  can  be  quite  large.  For  example,  for  planar  Ti5C4  and  Ti4C5 
clusters  consisting  of  9  atoms  the  number  of  geometrical  parameters  is  52  and  68, 
respectively.  To  optimize  structure  of  27-atom  cubic  cluster  of  TiuCn  930 
parameters  are  used. 
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Such  a  powerful  optimization  procedure  can  be  used  not  only  to  clarify  the 
configuration  of  atoms  in  the  clusters  already  complete,  but  also  in  the  analysis  of 
dynamic  behavior  of  free  atoms,  molecules,  and  small  fragments.  Some  of  these 
tasks  include  segmentation  of  several  nanoclusters,  separated  by  a  sufficiently 
large  distance  from  each  other.  This  section  examines  possible  cases  of 
segmentation  of  three  planar  nanoclusters  of  Ti5C4  and  Ti4C5  type,  spaced  at  some 
distance  in  two  coordinate  axes.  The  goal  of  computational  experiments  is  to 
address  two  issues.  First,  whether  segmentation  of  three  planes  into  cubic 
nanocluster  is  possible  and  under  what  conditions  can  this  happen?  Second,  how 
will  selected  planar  structures  behave  themselves  under  a  sufficiently  large 
disordering? 


- 

«  -^4 

v,f 

1 

/ 

tJ-Ai  ifl'* 

e  f  g  h 

a  -  initial  configuration;  b  -  step  25;  c  -  step  33;  d  -  step  36;  e  -  step  41;/-  step  42; 

g  -  step  50;  h  -  step  93,  optimized  structure 
Fig.  62.  Segmentation  of  planar  fragments  Ti5C4  -  T^Cs  -  Ti5C4  (Dy  =  0.4  nm,  Dz= 
0. 1  nm)  into  cubic  nanocluster  during  geometric  optimization 


In  a  simple  parallel  shift  of  three  planes  in  the  segmentation  into  cubic 
nanocluster  no  special  problems  do  emerge.  More  interesting  are  the  cases  of 


Distribution  A:  Approved  for  pubiic  reiease;  distribution  is  uniimited. 


mutual  movement  of  the  planes  along  two  axes.  The  results  of  five  segmentation 
variants  are  shown  in  Fig.  62-66.  Sequence  of  three  planar  fragments  Ti5C4  -  Ti4C5 
-  Ti5C4  we  denote  as  (1)  -  (2)  -  (3).  Origin  of  coordinates  is  placed  to  the  central 
atom  of  fragment  (2),  in  this  case  a  carbon  atom.  Fragment  (2)  lies  in  plane  XZ, 
fragments  (1)  and  (3)  are  separated  from  the  central  (2)  along  the  Y  axis  in  both 
directions.  Fig.  62  illustrates  an  optimization  process  of  the  first  variant  of  the 
system  with  fragments  shifted  along  the  Y  axis  by  an  Dy  =  0.4  nm  and  along  the  Z 
axis  by  Dz  =  0.1  nm. 

As  can  be  seen  from  Fig.  62  segmentation  of  planar  fragments  into  cubic 
nanocluster  Ti^Cia  (3x3x3)  is  a  rather  complex  process.  First  of  all,  we  note  that 
titanium  atoms  are  the  most  active,  it  is  the  metal-metal  interaction  of  adjacent 
planes  that  causes  the  attraction  of  fragments  (1)  and  (3)  to  the  center  (2).  After  the 
“skeleton”  is  built  with  titanium  atoms  (e),  a  process  of  ordering  of  carbon  atoms 
to  the  full  cube  symmetry  is  started.  Characteristically,  the  central  fragment  (2) 
remains  substantially  unchanged  and  the  picture  retains  symmetry  with  respect  to 
the  Y  axis,  because  planes  (1)  and  (3)  in  the  initial  configuration  are  shifted  with 
respect  to  (2)  by  the  same  distance. 

In  the  second  variant  the  segmentation  of  the  system  of  planar  fragments  with 
a  predominance  of  carbon  atoms  Ti4C5  -  Ti5C4  -  Ti4C5  into  cubic  nanocluster 
TiiaCu  (3x3x3)  was  considered.  The  coordinate  system  was  selected  similarly  to 
the  case  considered,  but  now  in  the  center  of  the  fragment  (2)  is  titanium  atom,  the 
displacements  along  the  axes  remain  the  same. 

The  optimization  process  proceeds  approximately  similar  to  the  scenario  of 
the  first  variant  but  with  some  difference.  At  the  first  stage  fragments  (1)  and  (3) 
approach  planar  fragment  (2)  and  attempt  to  establish  interaction  between  the 
metal  atoms.  Fig.  62c  and  Fig.  63c  illustrate  methods  of  bonding.  In  the  first  case, 
the  metal  atoms  of  the  central  chain  of  atoms  Ti  -  C  -  Ti  of  the  plane  (2)  interact 
with  the  metal  atoms  of  the  boundary  chains  Ti  -  C  -  Ti  of  fragments  (1)  and  (3) 
that  have  lost  their  flat  shape.  In  the  second,  on  the  contrary,  the  metal  atoms  of 
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central  chains  of  Ti  -  C  -  Ti  fragments  (1)  and  (3)  bind  to  the  metal  atoms  of  the 
boundary  chain  Ti  -  C  -  Ti  of  planar  fragment  (2).  But  this  interaction  is  unstable 
(Fig.  6M)  and  after  some  displacement  of  fragments  (1)  and  (3)  along  the  Z  axis  it 
becomes  possible  and  energetically  advantageous  to  establish  links  Ti  -  C  between 
the  boundaries  of  fragments  (Fig.  63e). 
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a  -  initial  configuration;  b  -  step  24;  c  -  step  25;  d  -  step  30;  e  -  step  40;/- 
step  70;  g  -  step  73;  h-  step  107,  optimized  structure 
Fig.  63.  Segmentation  of  planar  fragments  TuCs-  Ti5C4  -  T/Cs  (Dy  =  0.4  nm,  Dz 
=  0.1  nm)  in  a  cubic  nanocluster  during  geometrical  optimization 


A  similar  transition  was  observed  also  in  the  first  variant  (Fig.  62d,  e).  In  both 
cases,  these  bonds  were  formed  on  one  borders,  then  the  fragments  approached  and 
bonds  were  formed  on  the  opposite  boundaries  of  the  fragments  (Fig.  63/  g).  Then 
highly  symmetrical  cubic  shape  of  the  nanocluster  was  finally  established. 

In  the  third  variant  (Fig.  64)  the  segmentation  of  planar  fragments  Ti4C5  - 
Ti5C4  -  T/Cs  into  cubic  nanocluster  Ti^Cu  (3x3x3)  was  considered.  However, 
unlike  the  previous  case,  planes  (1)  and  (3)  are  shifted  with  respect  to  (2)  by  Dz  = 
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+0.1  nm  and  Dz  =  -  0.1  nm,  respectively.  Offset  along  X  axis  was  as  before  Dx  = 
0.4  nm. 


And  in  this  case  the  main  role  in  segmentation  is  played  by  the  interaction  of 
“metal  -  metal”  atoms  of  planes  (1)  and  (3)  with  atoms  of  (2).  The  central 
fragment  (2)  no  longer  maintains  the  shape  of  the  plane  due  to  displacement  of 
carbon  atoms  13  and  17  (Fig.  64,  b  and  c).  The  “skeleton”  of  cubic  nanocluster 
(Fig.  63,  d  and  e)  is  also  formed  mainly  by  long-range  bonds  of  titanium  atoms. 
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a  -  initial  configuration;  b  -  step  20;  c  -  step  30;  d  -  step  35;  e  -  step  45; 

/-  step  82  -  optimized  structure 

Fig.  64.  Segmentation  of  planar  fragments  TuCs  -  Ti5C4  -  T^Cs  (Dy  =  0.4  nm,  Dz 
=  0.1  nm)  in  a  cubic  nanocluster  during  geometric  optimization 


The  fourth  variant  with  the  system  Ti5C4  -  Ti4C5  -  Ti5C4  (Fig.  65)  differs 
from  the  first  in  that  the  central  plane  (2)  is  shifted  along  the  Z-axis  relative  to  the 
planes  (1)  and  (3)  by  0.2  nm.  Thus,  the  nearest  neighbors  of  the  upper  boundary  of 
planes  (1)  and  (3)  become  a  chain  of  atoms  Ti  -  C  -  Ti  in  the  center  of  plane  (2). 
And  in  the  first  steps  of  optimization  it  can  be  seen  (Fig.  65b,  c)  that  the  coupling 
of  fragments  (1)  and  (3)  takes  place  with  the  center  of  fragment  (2)  due  to  metal  - 
metal  interaction.  But  after  a  few  steps  atomic  configuration  changes  -  fragment 
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(2)  “is  pushed  out”  along  the  positive  direction  of  Z  axis  and,  as  can  be  seen  in  Fig. 
65g,  establishes  Ti  -  C  bonds  of  titanium  atoms  from  the  upper  boundaries  of 
fragments  (1)  and  (3)  with  carbon  atoms  of  the  lower  boundary  of  fragment  (2). 
Lower  boundaries  of  fragments  (1)  and  (3)  consist  of  the  same  type  of  atoms.  The 
attempt  to  establish  interaction  between  titanium  atoms  (Fig.  65e)  is  energetically 
less  advantageous  than  their  remote  location  (Fig.  65/,  g).  Optimized  shape 
resembles  the  bucket.  Central  fragment  (2)  retained  plane  shape  and  the  fragments 
(1)  and  (3)  arranged  symmetrically  with  respect  to  (2),  changed  their  shape 
significantly. 
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a  -  initial  configuration;  b  -  step  20;  c  -  step  25;  d  -  step  30;  e  -  step  37;/- 
step  45;  g  -  step  55\h-  step  77  -  optimized  structure 
Fig.  65.  Segmentation  of  planar  fragments  Ti5C4-  TuCs  -  Ti5C4  (Dy  =  0.4  nm,  Dz 
=  0.2  nm)  during  geometric  optimization 
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a  -  initial  configuration;  b  -  step  20;  c  -  step  25;  d  -  step  35;  e  -  step  40;/-  step  60; 

g  -  step  62\h-  step  85  -  optimized  structure 
Phc.  66.  Segmentation  of  planar  fragments  T^Cs  -  Ti5C4  -  TuCs  (Dx  =  0.4  nm,  Dy 
=  0.2  nm)  during  geometric  optimization 

Fifth  variant  of  segmentation  shown  in  Fig.  66,  has  an  initial  shape  of 
fragments  Ti4C5  -  Ti5C4  -  Ti4C5  similar  to  the  third  variant,  but  with  a  shift  of 
planes  (1)  and  (3)  with  respect  to  (2)  by  0.2  nm.  Shifting  planes  (1)  and  (3)  along 
axis  Z  by  a  distance  twice  as  much  has  led  to  the  configuration  where  the  nearest 
atoms  of  the  upper  boundary  (1)  and  the  middle  (2)  are  the  two  pairs  of  carbon 
atoms  and  only  one  -  of  titanium  atoms.  This  does  not  lead  to  the  formation  of 
bonds  between  them,  and  forces  (Fig.  66,  b)  fragments  (1)  and  (3)  move  toward  the 
upper  and  lower  boundaries  of  the  fragment  (2).  Initially,  metal  atoms  in  the 
middle  of  fragments  (1)  and  (3)  are  coupled  with  metal  atoms  at  upper  and  lower 
boundary  of  fragment  (2)  (Fig.  66c).  However,  this  configuration  is  unstable, 
established  bonds  are  broken  (Fig.  66d),  fragments  (1)  and  (3)  are  displaced  even 
more  (Fig.  66c)  and  a  more  advantageous  belt-type  configuration  is  formed.  The 
process  ends  up  by  rounding  of  free  boundaries  of  fragments  (1)  and  (3)  and  their 
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coupling  to  the  common  chains  of  fragments  (1)  -  (2)  and  (2)  -  (3).  The  last  stage 
of  the  segmentation  can  be  seen  in  Fig.  66/,  g,  h.  The  optimized  shape  resembles  a 
dimer,  it  is  symmetrical  relative  to  the  central  fragment  (2)  which  maintained  a 
planar  shape.  Fragments  (1)  and  (3)  rolled  up  into  distorted  rectangles  (2x1x1). 

The  conducted  computational  experiments  show  that  the  planar  Ti5C4  and 
Ti4C5  nanoclusters  under  certain  distance  conditions  can  spontaneously  segment 
into  cubic  nanoclusters  Ti^Cis  and  Ti^Cu.  The  process  of  assembly  takes  place 
exclusively  due  to  the  interatomic  interaction  of  9-atomic  fragments  with  each 
other.  Common  to  all  variants  is  a  long-range  metal  -  metal  interaction,  drawing 
fragments  to  each  other.  Only  after  the  fragments  have  shifted  so  much  that  can 
form  a  kind  of  “skeleton”,  the  process  begins  of  finding  the  most  advantageous 
from  the  energy  point  of  view  configuration  of  the  atoms.  In  the  case,  when  the 
distance  between  the  planes  is  about  0.4  nm,  which  is  close  to  twice  the  length  of 
the  Ti  -  C  bond  in  the  bulk  crystal,  and  a  mutual  shift  in  the  direction 
perpendicular  to  the  planes  is  in  the  range  of  ~  0. 1  nm  (about  half  the  length  of  the 
Ti  -  C  bond  in  the  bulk  crystal),  segmentation  of  planar  fragments  results  in  the 
formation  of  27-atomic  cubic  nanocluster.  With  a  larger  shift  (~  0.2  nm) 
segmentation  of  fragments  leads  to  non-cubic  forms. 

As  a  result  of  computational  experiments  a  number  of  possible  nanoclusters 
were  obtained:  two  nanoclusters  with  the  stoichiometric  composition  Ti^Co  with 
cubic  symmetry  (3x3x3)  (Fig.  62h)  and  in  the  form  of  a  “bucket”  (Fig.  65h),  as 
well  as  two  nanoclusters  Ti^Cu  with  the  structure  of  a  cube  (3x3x3)  (Fig.  63/ 
6Ah)  and  “dimer”  (Fig.  66h).  Table  20  shows  some  of  the  structural  and  energetic 
characteristics  of  these  nanoclusters.  Note  that  the  values  of  the  lengths  of  Ti-C 
bonds  in  all  four  variants  of  clusters  are  quite  close,  and  the  ratio  of  maximum  to 
minimum  bond  length  is  about  ~  1.08  -  1.09.  Only  for  “bucket”  cluster  is  this  ratio 
slightly  more  -  1.11.  The  disbalance  of  bonds  within  the  cluster  can  lead  to  the 
structural  instability.  In  this  sense,  the  least  stable  is  a  “bucket”  configuration. 
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Table  20.  Geometric  and  electron  properties  of  segmented  nanoclusters 


Characteristic 

Tii4Ci3 

“cube” 

TiuCn 

“cube” 

TiisCm 

“cube” 

TiibCm 

“AHMCp” 

Maximum  length  of  Ti-C  bond,  nm 

0.2158 

0.2088 

0.2194 

0.2101 

Minimum  length  of  Ti-C  bond,  nm 

0.1990 

0.1874 

0.2021 

0.1943 

Number  of  atoms  with  three  C-Ti  bonds 

0 

10 

8 

9 

Binding  energy,  eV/atom 

10.373 

9.510 

10.370 

9.804 

HOMO  -  LUMO,  eV 

0.666 

0.392 

0.716 

0.573 

Structural  stability  of  the  cluster  is  highly  dependent  on  the  amplitude  of 
oscillations  of  the  individual  atoms  or  group  of  atoms.  In  the  studied  objects  the 
most  mobile  are  light  carbon  atoms  that  have  three  or  four  bonds  with  titanium 
atoms.  Apparently,  it  is  the  triple  bonds  that  are  the  least  strong.  A  more  rigorous 
analysis  was  carried  out  in  [761  in  the  study  of  the  stability  of  cubic  nanoclusters 
Tii4Ci3,  found  experimentally  and  TioCu,  not  recorded  in  the  experiments.  The 
reason  of  instability  of  the  latter  is  that  the  cluster  contains  eight  carbon  atoms 
located  at  the  comers  of  a  cube  having  three  bonds  with  titanium  atoms.  In  this 
regard  it  is  useful  to  compare  for  the  clusters  considered  the  number  of  carbon 
atoms  having  three  C-Ti  bonds.  From  Table  20  it  can  be  seen  that  the  most  stable 
stmcture  is  cubic  nanocluster  TiuCo,  in  which  all  carbon  atoms  have  four  bonds 
with  the  atoms  of  titanium.  Cluster  “bucket”  has  8  carbon  atoms  with  triple  C-Ti 
bonds,  and  two  atoms  (angular  in  the  central  planar  fragment)  have  only  two  bonds 
with  atoms  of  titanium.  Estimates  of  the  magnitude  of  the  binding  energy  confirms 
the  findings:  highest  binding  force  is  in  cubic  structures  and  the  lowest  in  “bucket”. 

If  we  assume  that  considered  four  types  of  nanoclusters  are  stable  or 
metastable  and  may  participate  in  the  self-assembly  of  crystal  stmctures,  we  can 
point  out  the  possible  ways  of  further  segmentation  of  nanoclusters.  For  the  cubic 
nanoclusters  expansion  of  the  form  TiuCu  -  TioCu  in  all  three  directions  is 
obvious.  “Bucket”  nanocluster  can  be  stabilized  if  the  coupling  of  Ti^Cn  -  Ti^Cn 
will  go  on  in  such  a  way  that  the  planar  fragment  (2)  will  enter  into  the  gap 
between  the  fragments  (1)  and  (2)  and  this  will  lead  to  the  formation  of  new  Ti  -  C 
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bonds.  Planar  fragment  (2)  also  can  become  more  stable  if  contact  with  the 
fragments  (1)  and  (3)  will  take  place.  In  this  case,  double  C  -  Ti  bonds  will 
transform  to  triple.  Such  coupling  increases  the  deficit  of  carbon  atoms,  but 
crystalline  titanium  carbides  are  known  for  their  fairly  wide  range  of  homogeneity. 
Nanocluster  “dimer”  has  two  components  of  shape:  planar  fragment  (2)  and 
parallelepipeds  Ti6C6  (2x1x1)  formed  from  the  initial  fragments  (1)  and  (3). 

Ti6C6  nanoclusters  have  high  symmetry  and  can  be  stable  by  themselves. 
Calculations  of  binding  energy  and  of  atomic  oscillations  showed  that  IR  spectrum 
contains  no  imaginary  frequencies,  which  confirms  the  stability  of  the  Ti6C6 
cluster.  This  means  they  can  be  “building  blocks”  during  segmentation  into  the 
crystalline  structure.  Planar  fragment  (2)  has  probably  the  weakest  bonds,  where 
carbon  atoms  located  in  the  center  of  the  side  edges  each  have  three  C  -  Ti  bonds. 
Further  segmentation  and  stabilization  of  “dimers”  can  occur,  for  example,  when 
Ti6C6  nanoclusters  are  coupled  in  directions  that  allow  formation  of  new  C  -  Ti 
bonds.  Besides,  stabilization  of  fragment  (2)  is  possible  during  extension  of 
“dimer”  -  “dimer”  perpendicular  to  its  plane.  Such  coupling  leads  to  the 
enrichment  of  the  boundary  carbon  atoms  by  one  bond  with  titanium  atoms. 


Conclusion 

With  the  help  of  the  program  complex  “Gaussian'03”,  in  the  approximation 
of  density  functional  theory,  computational  experiments  on  the  segmentation  of 
carbide  nanoclusters  were  performed  for  the  first  time.  It  was  shown  that  even 
metastable  planar  forms  of  Ti5C4  and  Ti4C5  type  can  spontaneously  form 
nanoclusters  TiuCn  and  Tii^Cu  of  various  shapes  (cube,  “bucket”,  “dimer”), 
depending  on  the  distance  and  relative  orientation  of  free  surfaces.  Segmentation 
mechanism  is  due,  primarily,  to  long-range  interaction  of  metal  atoms  from 
different  planes.  After  the  formation  of  “skeleton”  from  metal  atoms,  formation  of 
Ti  -  C  bonds  and  final  configuration  takes  place. 
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III.  MEZOSCALE 


8.  Imitational  modeling  of  the  process  of  eutectic  co-crystallization  on  the 
basis  of  the  cellular  automaton  method. 

Eutectic  co-crystallization  of  boride-boride  composites  occurs  at  relatively 
high  (Tm  ~  2000-3000°C)  temperatures  and  is  associated  with  very  specific  physical 
conditions  [81].  This  necessitates  the  employment  of  the  formalism  of  nonequilibrium 
thermodynamics  and  nonlinear  equations  in  the  mathematical  models  developed  for 
the  description  of  these  processes.  In  the  nonlinear  thermodynamics,  there  is  a 
convenient  conventional  approach,  which  uses  the  partitioning  of  all  the  space  of  the 
subject  under  consideration  into  a  large  number  of  "elementary  volumes",  each  of 
them  being  specified  with  its  own  set  of  state  parameters.  Such  an  approach  is 
analogous  to  the  method  of  the  imitational  modeling  known  as  "cellular  automata" 
(CA)  [82].  At  present,  there  are  practically  no  works,  which  would  exploit  this 
methodology  for  simulating  real  processes  of  crystallization.  The  available 
publications  are  mostly  illustrative  and  aimed  at  the  demonstration  of  the  capabilities 
of  the  approach  rather  than  solving  actual  specific  problems  [83].  In  this  work  the  CA 
approach  is  implemented  for  2D  imitational  modeling  of  eutectic  structures. 

Developing  the  model 

In  order  to  define  the  "life"  rules  for  CA,  let  us  consider  the  simplified 
physical  model  of  the  crystallization  process  occurring  in  nonequilibrium 
dynamical  conditions. 

In  the  state  of  a  melt,  the  processes  of  nucleation  in  the  system  are  absent, 
as  are  the  conditions  needed  to  provide  the  crystal  growth.  The  melt  represents  the 
uniform  system  described  by  the  state  parameters,  which  determine  its  behavior. 
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In  the  process  of  cooling  the  system,  there  appear  local  deviations  of  the 
average  chemical  composition  (chemical  fluctuations),  or  local  changes  in  the 
temperature,  pressure,  etc.  (the  fluctuations  of  the  state  parameters  of  the  system). 

The  occurrence  of  the  local  meta- stability  in  the  system  gives  rise  to  the 
preconditions  for  the  nucleation  in  the  melt. 

Within  the  local  (physically  microscopic)  volume  in  the  meta-stable  state, 
i.e.,  far  from  equilibrium,  there  are  two  random  processes  possible,  namely,  the 
return  to  the  average  state  of  volume  within  the  system  ("dissolving"  the 
fluctuations  due  to  the  diffusion,  or  relaxation  of  fluctuations),  or  the  transition  to  a 
new  quality,  the  formation  of  a  solid  nucleus. 

The  crystal  growth  and  the  formation  of  the  inter-crystallite  boundaries 
occurs  from  the  "stable"  nuclei  of  the  leading  phase. 

On  this  basis,  we  can  formulate  the  following  stages  of  the  CA 
development  for  the  binary  system  to  be  modeled: 

-  onset  of  the  crystallization  process,  i.e.  the  appearance  of  a  fluctuation; 

-  diffusion  between  the  elementary  volumes; 

-  occurrence  of  the  initial  nuclei  of  the  eutectic  colony  in  the  two- 
dimensional  space; 

-  destruction  or  stabilization  of  nuclei; 

-  layered  crystal  growth; 

-  formation  of  inter-crystallite  boundaries; 

-  termination  of  the  process; 

-  formation  of  the  stable  (steady-state)  microstructure. 

To  perform  the  numerical  experiment,  we  chose  the  model  sample  based  on 
the  two-phase  LaB6-ZrB2  alloy  with  fiber  structure  (Figs.67,  68). 
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Fig.  67.  Schematic  of  the  composite 
to  be  modeled;  here  A  and  B  are  the 
components  mutually  insoluble  in  the 
solid  state,  A  is  the  matrix  (lanthanum 
hexaboride,  LaBe),  and  B  is  the 
fibered  diboride  phase  (zirconium 
diboride,  ZrB2). 


Fig.68.  Cross-section  microstructure  of 
the  directionally  crystallized  two- 
phase  alloy  LaB6-ZrB2  (micrograph 
obtained  by  scanning  electron 
microscope  «Stereoscan  S4-10»[84]. 


The  physical  parameters  of  the  sample  of  the  material  prototype  are  as 
follows:  [81]: 

•  eutectic  temperature  at  which  the  crystallization  occurs: 
Te  =  2740  +  40  K; 

•  eutectic  volume  fraction  for  ZrB2:  C  =  16%; 

•  average  fiber  diameter:  0,2  -  1,0  micron; 

•  inter-fiber  spacing:  0,8  -  1,2  micron. 

As  the  result  of  the  imitational  modeling,  the  CA  field  filled  with  the  cells  of 
two  types,  A  and  B,  should  be  formed,  with  the  relevant  structure  and  in 
accordance  with  the  above  physical  parameters.  Each  cell  type  has  to  be  specified 
via  the  vector  of  physical  and  chemical  material  parameters  (density,  specific  heat, 
latent  heat,  heat  and  temperature  conductivities,  etc.),  which  are  needed  for  further 
simulation. 
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The  simulation  has  to  be  performed  on  the  square  CA  field  consisting  of 
nxn  cells  (n  is  one  of  the  basic  model  parameters),  representing  the  cross-section 
of  the  sample.  Every  cell  may  occupy  one  of  three  states: 

•  empty  cell  (contains  0,  associated  with  the  melt); 

•  nucleus  cell  (contains  the  information  about  the  nucleus  and  fiber); 

•  fiber  cell  (contains  1,  associated  with  the  solid  phase  B). 

The  evolution  of  CA  field  proceeds  in  two  stages.  At  the  first  stage,  there 
emerge  B  phase  nuclei  on  the  CA  field.  At  the  second  stage,  they  grow  and  give 
rise  to  the  formation  of  fibers.  Within  the  proposed  model,  the  fiber  cross-section  is 
assumed  to  be  square. 

The  nucleus  cell  contains  all  the  information  related  to  the  nucleus  itself  and, 
partially,  to  the  fiber  originated  from  it.  The  fiber  cannot  exist  independently, 
without  a  nucleus.  The  nucleus  cell  contains  the  following  data: 

•  V,  y  —  nucleus  coordinates  in  the  field; 

•  status  —  nucleus  state;  can  take  three  values: 

o  “n”  — nucleus;  the  default  value, 
o  “g”  — growing  fiber, 
o  “f  ’  — full-grown  fiber; 

•  left,  right,  up,  down  —  the  number  of  cells,  which  the  nucleus  did  grow 
in  the  corresponding  direction. 

The  nucleus  can  grow  in  any  direction  to  form  the  fiber.  The  fiber  can  grow 
one  cell  left,  right,  up,  or  down  per  one  step,  while  retaining  the  rectangular  form. 
The  fiber  growth  is  limited  by  /x  /  size.  /  is  one  of  the  basic  model  parameters. 
The  examples  of  the  fiber  growth  are  given  in  Fig.  69.  Alternatively,  the  nucleus 
can  "die".  The  fiber  grown  from  it  "dies"  therewith  as  well;  all  the  related 
information  is  deleted  and  all  fiber  cells  are  filled  with  zeros. 
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Fig.  69.  Two  possible  implementations  of  the  successive  states  in  the  process  of 
fiber  growth  up  to  the  3x3  cell  size.  The  nucleus  and  the  growing  fiber  are  shown 
in  black  and  grey,  respectively. 

Some  neighborhood  near  every  nucleus  or  fiber  is  viewed  as  the  "dead 
zone".  The  dead  zone  is  defined  by  the  following  considerations:  a)  within  the  dead 
zone,  no  other  nucleus  can  emerge;  b)  growing  fibers  cannot  occupy  the  dead  zone 
cells  belonging  the  other  fibers;  c)  dead  zone  width  g  is  one  of  the  basic  model 
parameters.  Within  the  framework  of  the  model,  it  is  allowed  to  set  different  dead 
zone  widths  on  the  stages  of  nucleation  ( )  and  fiber  growth  ( ). 

A  nucleus  neighborhood  of  the  width  g  is  defined  as  the  set  of  cells  with  the 
coordinates,  which  differ  from  the  nucleus  coordinates  no  more  than  by  g  in  both 
directions.  Therefore,  the  neighborhood  of  the  nucleus  with  the  coordinates 
contains  all  cells  with  the  coordinates  (x,  j)  satisfying  the  relation 

Xo-g^X<X„  +  g 

yo-g^y^yo  +  g 

The  fiber  neighborhood  of  a  radius  g  includes  all  cells  with  the  coordinates, 
which  differ  from  the  coordinates  of  the  nearest  fiber  cell  no  more  than  by  g . 
Therefore,  for  the  case  that  the  nucleus  has  the  coordinates  and  did  grow  by 
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left,  right,  up,  down  cells  in  each  of  these  directions,  the  neighborhood  will  be 

represented  by  the  cells  with  the  coordinates  (x,  y)  fulfilling  the  relation 

X,)  -  left  -  ^  <  X  <  Xo  +  right  +  g 
yQ-up-g<y<y^+down  +  g 

An  example  of  the  nucleus  and  fiber  neighborhood  is  given  in  Fig.  70.  In 
this  model,  the  notion  of  the  neighborhood  is  introduced  for  nucleus  and  fiber  only. 
The  neighborhood  of  an  arbitrary  cell  is  undefined. 


a)  b) 

Fig.  70.  Example  of  the  nucleus  (a)  and  fiber  (b)  neighborhood.  The  nucleus  and 
the  fiber  are  displayed  in  black  and  their  neighborhood  in  grey,  for  =  g^=?>. 

Initially,  the  CA  field  is  filled  with  zeros. 

At  the  first  stage,  the  nucleation  on  the  CA  field  is  modeled;  arf  attempts  of 
seeding  are  made  (a  is  one  of  the  basic  parameters  of  the  model).  The  first  stage  of 
nucleation  proceeds  as  follows:  a  random  cell  from  the  set,  where  the  seeding  yet 
not  been  attempted,  is  chosen;  then,  in  this  cell  a  nucleus  is  created.  In  the  case  that 
the  nucleus  appears  in  the  dead  zone  of  another  nucleus,  it  dies  immediately. 

At  the  second  stage,  the  revision  of  all  the  nuclei  in  the  order  of  seeding  is 
made.  Such  a  revision  repeats  2/  -2  times.  In  the  course  of  revision,  an  attempt  for 
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every  nucleus  to  grow  by  one  cell  is  made.  The  direction  of  growth  is  chosen  at 
random,  with  the  following  rules  being  fulfilled: 

•  as  the  result  of  the  growth,  the  fiber  should  not  occupy  the  dead  zone  of 
another  fiber  or  nucleus; 

•  it  should  not  leave  the  area  of  CA  field; 

•  its  size  in  both  directions  cannot  exceed/  cells. 

In  the  case  that  there  are  no  directions  allowed  for  fiber  growth,  which 
would  meet  above  conditions,  and,  at  the  same  time,  the  fiber  did  not  reach  its 
maximum /x/  size,  it  dies  and  is  excluded  from  further  revisions. 

As  the  result  of  this  process,  there  forms  a  number  of  square  fibers  of 
/x/cell  size  on  the  CA  field.  By  using  the  number  of  the  grown  fibers,  one  can 


evaluate  the  concentration  for  the  phase  B  as  the  ratio  of  the  number  of  cells 
occupied  by  fibers  to  the  total  number  of  cells. 


(83) 


Description  of  the  program 

To  implement  the  model,  a  program  has  been  written  by  using  the  Python 
programming  language.  Pseudo-random  number  were  generated  by  using  the  built- 
in  generator  "Mersenne  twister".  To  visualize  the  cellular  automaton,  the  well 
known  library  Matplotlib  has  been  exploited. 

Basic  parameters  of  the  model  and  program  are  as  follows: 

•  n  —  field  size. 

•  a  —  the  ratio  of  the  number  of  seeding  attempts  to  the  total  cell  number; 
in  total,  an  seeding  attempts  are  made. 

•  /  —  size  of  the  full-grown  fiber. 

•  g  —  dead  zone  width: 
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O 


in  the  process  of  fiber  formation; 


o  —  in  the  process  of  fiber  growth. 

The  program  is  implemented  as  the  library  written  in  Python  language.  It  has 
command  line  (not  graphical)  interface.  The  model  can  be  launched  by  the  regular 
library  call,  for  instance, 

import  sys 

sys. path. append  {"[path  to  the  "models”  directory]") 

from  modelb  import  * 

Model6(n  =  300,  a  =  0.1,  f  =  8,  gl  =  14,  g2  =  7).run(step  =  0.1,  plot  =  True) 

Such  an  approach  simplifies  the  numerical  experiment,  because  it  enables 
one  to  start  the  model  with  different  parameters  and  accumulate  the  results  of  the 
work.  Alternatively,  the  program  can  be  started  by  directly  running  the  source  code 
file  modelb.py;  in  this  case  the  parameters  have  to  be  specified  within  the  file  itself 
(Fig.  71). 

The  process  of  nucleation  can  take  long  time.  To  indicate  the  progress,  the 
program  displays  the  intermediate  results  (the  number  of  attempts  and  the  nuclei 
seeded)  with  the  period  determined  by  the  parameter  'step':  the  period  is  equal  to 
atd  ■  step .  At  the  stage  of  fiber  growth,  the  program  displays  the  intermediate  results 
after  each  revision,  namely,  the  current  revision  number  and  the  current  number  of 
the  fibers  left.  Then,  the  concentration  of  the  fiber  material  is  calculated  in 
accordance  with  Eq.  (83)  and  the  field  image  is  plotted.  An  example  of  the 
program  output  is  given  in  Fig.  7 1 . 
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Fig.  71.  The  example  of  the  program  output  started  from  the  source  code  editor 
Sublime  Text. 

The  program  user  has  to  specify,  in  first  place,  the  value  of  the  fiber 
concentration.  For  the  eutectic  composition  of  LaB6-ZrB2  ,  the  latter  amounts  to 
16%.  After  that,  the  initial  values  for  basic  parameters  are  set  and  the  program  gets 
launched.  A  good  choice  for  the  initial  values  are  the  ones  given  in  Fig.  7 1 .  If  the 
resulting  fiber  concentration  is  far  from  the  desired  value,  the  basic  parameters 
have  to  be  adjusted  until  the  needed  concentration  is  obtained.  First  of  all,  the 
parameters  / ,  and  should  be  varied.  As  was  observed  in  practice,  the  value 


for  g^  must  be  close  to  g^-f  +  l.  Otherwise,  the  image  of  CA  field  gets  essentially 


different  from  that  given  in  Figs.  67  and  68. 
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The  program  has  the  option  to  vary  the  field  size  parameter  n .  With 
increasing  n ,  the  computational  time  grows  rapidly.  In  order  to  reduce  it,  one  can 
diminish  the  parameter  a,  which  affects  the  number  of  seeding  attempts.  In  so 
doing,  it  is  desirable  to  keep  the  number  of  nucleation  attempts  an  much  greater 
than  the  maximum  number  of  nuclei  to  be  seeded  (the  latter  obtains  at  a  =  1.0). 
Otherwise,  the  structure  of  CA  field  obtained  in  simulation  will  be  different  from 
the  one  needed.  By  using  the  personal  computer  of  typical  configuration,  the 
reasonable  computation  time  of  the  program  corresponds  to  the  numbers  of  n 
within  the  limit  of  1000  cells. 

Testing  the  program  by  varying  the  basic  parameters 

For  testing  purposes,  the  numerical  experiment  has  need  performed. 

In  particular,  the  effects  of  basic  parameters  on  the  final  B  phase 
concentration,  the  uniformity  of  fiber  distribution,  and  the  computational  time  has 
been  examined. 

Study  of  the  boundary  effects 

According  to  the  rules  accepted  within  this  model,  there  exists  the  minimum 
distance  between  nuclei/fibers,  which  is  equal  to  g  cells.  At  the  same  time,  the 
rules  do  not  limit  the  minimum  distance  between  nuclei/fibers  and  the  field 
boundary.  This  suggests  that  the  fiber  concentration  near  field  boundaries  can  be 
higher  than  that  at  distance.  This  assumption  is  confirmed  by  the  results  of 
numerical  experiment  with  varying  the  field  size  (see  below).  The  elevated  fiber 
concentration  near  field  boundary  is  demonstrated  by  Fig.  7 1 . 

In  the  numerical  experiment,  the  following  notations  were  used  (Fig.  72): 

•  d  — distance  from  the  given  cell  to  the  field  boundary,  exactly,  the 
minimum  distance  between  the  cell  and  field  boundaries  expressed  in  the 
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number  of  cells.  Therefore,  for  the  cells  located  immediately  on  the  field 
boundary,  J  =  0  (cells). 

•  Layer  of  cells  at  the  separation  d  from  the  boundary  is  the  set  of  all 
cells,  for  which  d  has  the  same  value  (i.e.,  equally  separated  from  the 
field  boundary). 

The  fiber  concentration  in  the  given  layer  is  defined  as  the  ratio  of  the 
number  of  cells  in  this  layer  occupied  by  the  fibers  to  the  total  number  of  layer 
cells.  The  program  has  been  restarted  100  times  with  the  following  parameters: 
n  =600,  a  =0.04,  /  =14,  .gi=30,  =15.  In  every  run,  a  new  field  was  obtained.  For 

each  field,  the  fiber  concentration  in  each  layer  has  been  evaluated.  After  that,  the 
concentration  C  averaged  over  all  the  results  for  the  given  d  was  calculated.  The 
dependence  of  the  average  concentration  C  on  the  distance  d  is  given  in  Fig.  73. 
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a)  b) 

Fig.  72.  Cell  layers  located  at  the  given  separation  from  the  field  boundary,  by 
taking  as  an  example  the  10  x  10  cell  field,  a)  the  layers  distinguish  by  colors,  the 
values  for  d  are  given  in  cells,  b)  colors  indicates  the  possible  fiber  configuration. 
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c,% 


Fig.  73.  Dependence  of  the  B  phase  concentration  on  the  separation  from  the  field 
boundary. 

The  numerical  experiment  has  confirmed  the  existence  of  the  boundary 
effect.  Considerable  increase  in  concentration  is  observed  near  the  field  boundary. 
Starting  from  the  field  boundary,  the  concentration  grows  with  the  parameter  d 
and  reaches  its  maximum,  which  is  1,5  -  2  times  greater  than  the  average  value. 
The  peak  is  located  approximately  at  the  distance  /  from  the  field  boundary.  The 
peak  is  followed  by  the  "gap",  i.e.  the  fiber  concentration  reduces  and  becomes  less 
than  the  average  value,  however,  with  the  deviation  much  less  than  that  of  the 
peak.  Further  oscillations  do  not  exceed  the  standard  deviation  for  the 
concentration.  The  standard  deviation  increases  while  approaching  the  field  center. 
Apparently,  this  is  due  to  the  reduce  in  the  number  of  cells  in  the  relevant  layers. 
Thus,  the  layers  which  are  close  to  the  field  center  should  be  dropped  in  analysis. 

Similar  results  were  obtained  for  the  other  parameters  (n  =200,  a  =1.0,  /  =4, 
^1=7,  ^2-3,  50  runs).  The  ways  to  diminish  or  eliminate  the  boundary  effects 
require  further  refinements  of  the  model  proposed. 
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Investigation  of  the  concentration  dependence  on  the  field  size 


To  study  the  dependence  of  the  B  phase  concentration  on  the  field  size  n ,  a 
series  of  numerical  experiments  has  been  performed.  The  parameters  were:  a  =  02, 
/  =  8,  =  14 ,  g^=l .  The  field  size  n  has  been  varied  within  the  range  from  50  up 

to  1000  cells.  For  each  value  of  n ,  there  were  made  10  program  runs,  and  the 
results  were  averaged.  The  dependence  obtained  is  given  in  Fig.  74. 


Fig.  74.  Dependence  of  the  B  phase  concentration  on  the  field  size. 

The  results  of  the  numerical  experiment  have  shown  that  the  increase  in  the 
field  size  results  in  the  moderate  reduce  of  the  B  phase  concentration.  For  larger 
field  sizes,  B  phase  concentration  is  nearly  constant;  the  changes  do  not  exceed  the 
standard  deviation.  The  standard  deviation  of  the  concentration  decreases  rapidly 
with  growing  n  (presumably,  in  accordance  with  the  hyberbolic  law). 

Study  of  the  dependence  of  the  computational  time  of  the  program  on  the  field  size 

To  examine  the  temporal  computational  complexity  of  the  algorithm,  a 
series  of  numerical  experiments  with  the  parameters  a  =1,0,  /=  8,  .gi=15,  g^=l 
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was  performed.  The  field  size  n  was  varied  within  the  range  from  50  to  800  cells 
with  adjusting  the  values  of  n  in  such  a  way,  that  every  next  n  was  approximately 
^  times  greater  than  the  foregoing  one.  For  every  n  within  the  limits  from  50  to 
400,  10  runs  were  performed,  and  for  the  values  500,  630  and  800  the  program 
was  run  only  once  (per  one  value)  in  view  of  considerable  increase  of 
computational  time. 

The  operation  time  of  the  program  was  determined  by  using  the  four  core 
CPU  AMD  Athlon^^  II  X4  620  at  the  frequency  2,60  MHz  (with  the  use  of  one 
core  only).  The  averaged  results  of  measurements  are  given  in  Table  21  and  in 
Fig.75.  The  analytical  approximation  for  the  results  of  measurements  has  been 
derived,  which  indicates  that  the  obtained  dependence  can  be  fitted  well  by  the 
following  relation: 

t  =  (84) 

where  A  =0,2487,  B  =-0,5509,  C  =  9, 1825  10“^  c. 

As  is  seen  from  Fig.  75  and  Table  21,  Eq.  (84)  approximates  the 
experimental  data  rather  well. 


a)  b) 


Fig.  75.  Dependence  of  the  operation  time  of  the  program  on  the  field  size  (the 
experimental  data  and  their  approximation  according  to  Eq.(84))  on  linear  (a)  and 
logarithmic  (b)  scale. 
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Table  21 


a)  comparison  of  the  measured  time  with  the  fit; 

b)  predictions  for  the  time  by  using  Eq.  (84)  for  the  numbers  in  the  range 
from  1000  to  10000. 


n,  cells 

prediction 
for  t,  s 

prediction  for  t 
in  alternative 
units 

1000 

5568,428 

93  minutes 

2000 

140268,9 

39  hours 

3000 

1086524 

12,6  days 

4000 

4987953 

2  months 

5000 

16946530 

6,5  months 

6000 

47270274 

1,5  year 

7000 

1,15-10* 

3,6  years 

8000 

2,50-10* 

8  years 

9000 

5,04-10* 

16  years 

10000 

9,51-10* 

30  years 

n,  cells 

t,  s 

approximation 
for  t,  s 

50 

0,26 

0,258059 

63 

0,44 

0,443112 

80 

0,82 

0,806692 

100 

1,45 

1,464924 

125 

2,73 

2,757022 

160 

5,76 

5,785399 

200 

11,87 

11,73934 

250 

24,56 

24,68721 

320 

59,03 

58,57217 

400 

133,36 

132,8021 

500 

318,7 

312,0586 

630 

USA 

785,0108 

800 

2113,3 

2120,792 

a)  b) 


The  derived  relation  (84)  indicates  that  the  operation  time  grows  drastically 
with  increasing  the  field  size.  For  the  field  size  as  much  as  several  thousand  cells, 
the  operation  time  becomes  impractical.  The  increase  in  computing  capacity  cannot 
remedy  the  problem.  For  instance,  for  n>  1000  cells,  and  for  the  given  operation 
time,  the  1000-fold  increase  in  computing  capacity  enables  one  to  increase  the  field 
size  by  less  than  one  order.  One  should  expect  that  changing  the  parameters  will 
yield  the  similar  results.  Solution  to  the  problem  could  be,  probably,  the 
optimization  of  the  algorithm  of  nucleation. 
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Study  of  the  dependence  of  operation  time  of  the  program  on  the  parameter  a. 

As  was  demonstrated  by  the  above  mentioned  experiment,  the  increase  in 
field  size  results  in  rapid  increase  in  simulation  time.  The  examination  of  the 
functioning  of  the  program  has  revealed  that  the  most  computational  efforts  are 
spent  for  simulations  of  the  process  of  nucleation,  whereas  the  nucleus  growth  is 
much  less  time  consuming  (the  latter  proceeds  2-3  orders  more  rapidly).  Therefore, 
it  is  the  nucleation  process  that  should  be  optimized  to  make  simulations  more 
efficient.  For  instance,  the  nuclei  can  be  generated  in  randomly  selected  cells  only, 
instead  of  seeding  every  cell.  For  this  purpose,  we  introduced  the  parameter 
a  (0<  a  <  1)  defined  as  the  fraction  of  the  cells,  wherein  the  seeding  was  attempted. 

Thus,  the  nucleation  requires  an^  attempts  to  be  made.  To  make  the  simulation 
time  reasonable,  it  is  advisable  to  reduce  the  value  of  a  with  increasing  the  field 
size.  It  is  natural  to  expect  that  the  simulation  time  depends  linearly  on  a 
(provided  that  the  other  parameters  remain  unchanged).  To  verify  this  assumption, 
the  numerical  experiment  has  been  performed.  The  program  was  started  with  the 
following  parameters:  n  =  250,  /  =  8,  g^  =  15,  =1.  The  parameter  a  varied 

from  0.0  till  1.0  with  the  increment  of  0.1.  For  each  a,  the  program  has  been 
repeatedly  started  20  times.  The  averaged  experimental  results  were  used  to  plot 
the  relevant  dependence  along  with  the  approximation  by  the  linear  function 
(Fig.76).  In  doing  so,  the  point  ?(0)=0,04994c  was  regarded  as  being  known  and 

was  measured  independently. 
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Fig.  76.  Dependence  of  the  operation  time  of  the  program  on  the  parameter  a  along 
with  the  linear  approximation. 

As  is  seen  from  Fig.  76,  the  dependence  t[a)  for  given  model  parameters  is 

well  reproduced  by  the  linear  function  starting  from  the  point  near  the  coordinate 
origin.  It  is  natural  to  expect  that  the  other  parameters  will  yield  similar  results. 
Therefore,  the  model  operation  time  can  be  assumed  to  be  proportional  to  the 
parameter  a . 

It  should  be  noted,  that  diminishing  the  parameter  a  results  in  decrease  in 
the  number  of  nuclei,  and,  as  the  consequence,  in  decrease  of  B  phase 
concentration.  To  form  the  field  with  given  concentration,  the  other 
parameters (/,  ^1,^2)  should  be  adjusted. 

Distribution  of  the  number  of  nuclei 

The  results  of  numerical  experiments  for  the  same  set  of  parameters  differ 
because  of  the  use  of  the  pseudo-random  number  generator  in  simulations.  To 
examine  the  distribution  of  the  number  of  nuclei  (at  the  nucleation  stage;  the 
growth  stage  was  not  considered),  the  following  numerical  experiment  was 
performed.  The  program  has  been  restarted  repeatedly  with  the  parameters  n  = 
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100,  a  =  1,0,  /  =  4,  =  7,  and  the  number  of  nuclei  created  was  monitored.  The 

total  number  of  runs  was  3238,  and  the  results  were  used  to  plot  the  relevant  bar 
chart  (Fig.77).  The  following  quantities  were  evaluated: 

•  Average  number  of  nuclei:  103,5 1 ; 

•  Standard  deviation:  3,35; 

•  Coefficient  of  asymmetry:  0,0343; 

•  Coefficient  of  excess:  -0,0693. 

The  curve  for  the  normal  distribution  (scaled  by  the  factor  3238)  with  the 
same  parameters  is  plotted  in  Fig.  77  as  well.  Fig.  77  demonstrates  that  the 
distribution  of  the  number  of  nuclei  is  close  to  the  normal  one. 


Fig.  77.  Bar  chart  for  the  distribution  of  the  number  of  the  nuclei  seeded  in  3238 
runs.  For  comparison,  the  Gaussian  distribution  (scaled  by  3238)  is  displayed. 


Visualization  of  the  results  of  the  numerical  experiment 

2D-visualization 

In  the  present  model,  the  fiber  cross-section  is  assumed  to  be  square.  For  the 
purpose  of  visualization,  the  question  was  considered  as  of  replacing  the  squares 
by  the  circles  of  equal  area,  with  the  coinciding  centers. 


Distribution  A:  Approved  for  pubiic  reiease;  distribution  is  uniimited. 


This  feature  was  implemented  by  the  Matplotlib  library.  An  example  for  the 
program  output  for  the  parameters  n  =  300,  <2  =  0,1,  f=  8,  gi  =  14,  g2  =  l  is  shown 
in  Figs.  78,  79.  In  Fig.  78,  the  CA  field  with  the  square  fibers  is  displayed.  In  Fig. 
79,  the  same  field  is  represented  by  the  fibers  with  circle  cross-sections.  For  better 
visualization  of  circles,  the  field  size  was  scaled  from  300  x  300  to  1500  x  1500, 
i.e.  each  single  cell  of  the  initial  field  was  replaced  by  25  (5  x  5)  cells.  A 
possibility  is  provided  to  specify  the  different  field  size,  by  using  the  optional 
scaling  parameter  z  (which,  in  the  case  under  consideration,  is  equal  to  5). 


0  50  100  150  200  250  0  200  400  600  800  1000  1200  1400 


Fig.  78.  Initial  CA  field.  Fiber  cross-  Pnc.  79.  Initial  CA  field.  Fiber  cross- 
section  is  square.  section  is  circle. 

3D  visualization 

For  3D  visualization,  we  employed  the  Mayavi  application,  which  can  be 
used  as  a  library  and  as  an  independent  application. 

By  using  Mayavi,  starting  from  the  initial  CA  field  with  the  circle  cross- 
sections  of  fibers  (Fig.  79),  the  cylinders  can  be  constructed,  which  imitate  the 
configuration  of  fibers  in  a  real  sample  (Figs.  80,  81).  The  height  of  the  sample  can 
be  adjusted  by  the  parameter  h. 
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Phc.  80.  3D  visualization  of  the  Pnc.  81.  3D  visualization  of  the  initial 
initial  CA  field.  CA  field  with  plotting  sample  contour. 


The  Mayavi  application  has  the  option  of  changing  the  orientation  of  the  3D 
model  (Fig.  82). 


Phc.  82.  3D  visualization  of  the  initial  CA  field  with  plotting  sample  contour. 
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Conclusions 

A  model  for  imitational  2D  simulation  of  eutectic  structure  formation  based  on 
the  CA  technique  has  been  developed  and  implemented  as  a  computer  program.  As 
the  model  sample,  the  two-phase  alloy  LaB6-ZrB2  was  employed.  The  program 
implementation  of  the  proposed  model  is  described.  The  basic  parameters  of  the 
model  are  determined,  and  their  effects  on  the  structure  formation  are  investigated.  By 
means  of  numerical  experiments,  the  ranges  for  the  basic  model  parameters,  which 
provide  the  adequate  imitation  of  the  eutectic  structure  formation,  are  established. 


9.  Fractal  characteristics  of  the  destruction  surfaces 

In  the  investigation  of  the  deformation  and  destruction  of  materials,  one  of 
the  important  characteristics  is  the  structure  of  the  material  itself.  To  describe  the 
latter,  in  recent  years,  the  formalism  of  the  fractal  geometry  is  employed  [85-871. 

However,  at  present  there  are  no  unified  algorithms  for  calculations  of 
fractal  characteristics  of  materials  available,  in  particular,  for  images  obtained  in 
various  ways,  such  as,  for  instance,  via  the  electron  microscopy. 

In  our  earlier  investigations  [881,  we  have  shown  that  the  programs  widely 
used  for  evaluation  of  the  fractal  dimensionality  have  a  number  of  drawbacks. 
These  are,  in  first  place,  the  high  instability  and  sensitivity  to  image  size. 

In  view  of  this,  we  have  implemented  the  algorithm  for  calculation  of  the 
block  (Hausdorff)  fractal  dimensionality  and  performed  the  investigation  of  the 
behavior  of  the  program  while  processing  both  the  artificially  designed  model 
images  with  pre-defined  characteristics  like  Sierpinski  carpet,  and  real  images 
obtained  by  means  of  electron  microscopy. 

On  the  basis  of  the  approach  employed,  an  algorithm  for  calculations  of  the 
multifractal  characteristics  of  two-dimensional  images  has  been  developed  and 
implemented. 
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Algorithm  for  calculation  of  the  fractal  dimensionality  of  two-dimensional 
images 

The  analysis  based  on  the  formalism  of  fractal  geometry  plays  an 
increasingly  important  role  among  the  mathematical  models  [85-87,  89-901. 

Since  the  1980s,  the  notion  of  fractal  structures,  or  fractals,  gained  wide 
acceptance  among  researchers. 

The  entities,  which  can  be  characterized  by  fractional  dimensionality  of  their 
spatial  components,  are  called  fractals.  The  objects  with  the  fractional  geometrical 
dimensionality  gained  the  attention  of  mathematicians  since  the  beginning  of  20th 
century,  however,  they  did  not  find  any  practical  application  in  physical  and 
chemical  investigations.  The  interest  in  fractal  structures  has  been  awakened  after 
the  publication  of  Mandelbrot's  monographs  [91,  92],  where  it  was  shown  that  the 
vast  majority  of  natural  systems  can  be  regarded  as  fractals.  The  adoption  of  the 
notion  of  fractals  made  possible  a  deeper  insight  into  the  structural  organization  of 
various  physical  and  chemical  systems,  as  well  as  the  regularities  of  diverse 
processes  occurring  therein. 

There  are  several  ways  of  defining  the  fractal  dimensionality.  From  the 
viewpoint  of  mathematics,  these  definitions  are  not  equivalent.  The  details 
concerning  the  introduction  of  the  fractal  dimensionality,  can  be  found  in  one  of 
the  most  cited  monographs  in  this  area  [87]. 

The  evaluation  of  the  block  fractal  dimensionality  can  be  described  briefly 
as  follows.  Consider  a  set  P  of  points  in  Euclidean  space  with  the  dimensionality 
E  (in  the  case  of  plane  images  E  =  2).  The  fractal  dimensionality  of  P  is 
evaluated  by  means  of  the  following  procedure.  Let  us  cover  the  set  of  points  by 
the  net  of  E-  dimensional  cells  of  the  size  Ar ,  and  denote  the  number  of  the  cells 
containing  at  least  one  element  of  the  set  P  as  A/^(Ar).  The  sought-for  block 
fractal  dimensionality  of  the  set  can  be  derived  from  the  relation 
iV(Ar)D 
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This  approach  for  defining  the  fractal  dimensionality  directly  relates  to  the 
classical  definition  of  fractal  dimensionality,  which  is  accepted  as  the  basic  one  in 
[93]:  the  fractal  dimensionality  of  any  physical  object  can  be  found  from  the 


InA^(^) 

relation  =lim — ,  ,  ,  , 


where  A^is  the  number  of  blocks  (spheres)  of  typical 


linear  size  6,  which  completely  cover  the  subject  under  consideration. 

In  the  case  of  a  two-dimensional  black  and  white  image  (which  is  the  case 
typical  for  electron  microscopy),  the  set  P  is  represented  by  the  image  itself,  and 
the  net  of  cells  is  nothing  but  the  partitioning  of  the  image  into  the  squares.  Since 
the  simplest  way  to  obtain  the  fractal  dimensionality  is  the  least  square 

technique  applicable  to  the  linearized  expression  \nNU  -Z)5ln(Ar),  it  is  desirable 
that  the  values  of  In(Av)  would  be  distributed  uniformly.  In  other  words,  if  we  are 
going  to  obtain  M  experimental  points  of  the  function  In  =  /  (In  Ax) ,  then  the 
relevant  values  for  Ax  should  be 


Ax,.  = 


jn^. 

,M-1 


i=0,M-l, 


(85) 


where  S  is  the  maximum  image  size  in  pixels,  M  is  the  total  number  of  points. 
Notice  that  the  actual  total  number  of  points  may  be  less  than  M  in  the  case  that 
rounding  off  will  yield  equal  values  of  Ax,  for  different  i. 

Technically,  it  is  more  efficient  to  consider  the  quantity  2-Dj^,  i.e.,  in  order 


to  determine  the  fractal  dimensionality,  the  dependence  A/^(Ax)Ax^  □  Ax^  (or,  in 


linearized  form.  In NAx^  □  (2-Z)g)lnAx)  should  be  employed. 
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Computer  implementation  of  the  algorithm  and  the  study  of  its  properties 

Description  of  the  program  for  determining  the  block  fractal  dimensionality  of 
two-dimensional  images 

The  program  for  processing  the  images  obtained  from  the  electron 
microscope  based  on  the  above  technique  for  determining  the  fractal 
dimensionality  is  written  in  the  programming  language  Open  Watcom  C++ 
(http://www.openwatcom.org)  with  the  use  of  the  library  wxWidgets 
(http://www.wxwidgets.org).  Both  of  these  products  are  freely  distributed  as  open- 
source  programs;  in  addition,  they  both  allow  the  compilation  of  the  applications 
for  Windows,  OS/2  and  Linux  operating  systems.  Thus,  the  program  developed 
can  be,  in  principle,  transferred  to  any  of  the  above  listed  operating  systems.  At 
present,  it  operates  under  the  system  MS  Windows  (it  is  verified  and  tested  under 
MS  Windows  XP  n  Windows  7,  although  there  are  no  any  specific  requirements 
for  operating  systems). 

Presently,  the  developed  program  provides  the  evaluation  of  both  block  and 
mass  fractal  dimensionality.  The  overall  appearance  of  the  program  window  is 
given  in  Fig.  83. 

The  program  enables  one  to  set,  prior  to  computations,  the  global  tunings 
(which  are  saved  in  registry  file  and,  therefore,  remain  valid  for  subsequent 
program  runs),  as  is  shown  in  Fig.  84. 

The  first  global  setting  is  the  selection  of  the  algorithm  (block  or  mass)  for 
calculations,  the  second  one  is  the  number  of  points  of  experimental  curve.  The 
calculation  of  the  derivative  of  the  experimental  curve  (which  is,  in  fact,  the  basic 
issue  for  the  evaluation  of  fractal  dimensionality)  can  be  based  on  the  algorithms 
for  calculation  of  the  derivative  by  means  of  constructing  the  finite-difference 
scheme,  with  the  use  of  three  or  five  points.  The  last  parameter  specifies  the  units, 
which  are  used  to  set  the  threshold  for  transforming  the  images  into  the  black  and 
white  mode  (the  meaning  of  this  parameter  will  be  considered  in  details  below). 
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Fig.83.  Overall  look  of  the  window  of  the  program  for  computing  the  fractal 
dimensionality  of  images. 
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Fig.84.  Window  of  global  settings  of  the  program  for  computing  the  fractal 
dimensionality  of  images. 
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The  program  works  as  follows.  After  selecting  the  global  settings  (or,  if  you  don't 
open  the  settings  window,  accepting  the  previous  ones),  you  open  the  image  file 
(the  basic  graphical  formats,  BMP,  GIF,  JPEG  and  TIFF,  are  supported).  If  the  file 
relates  to  the  color  image,  it  transforms  to  the  grayscale  format  according  to  the 
conventional  algorithm  implying  that  a  pixel  with  red,  green,  and  blue, 
accordingly,  R,  G,  and  B,  is  replaced  by  the  pixel  from  grayscale  with  the  value 
0.299xR+0.587xG+0.114xfi  (as  recommended  by  ITU-R  BT.601).  A  bar  chart 
indicating  the  number  of  pixels  with  different  brightness  is  displayed. 

Then  the  image  is  transformed  to  black  and  white  mode  (1  bit  per  pixel), 
which  provides  the  basis  for  further  computations.  Here  one  can  control  the 
process  of  transformation  by  setting  the  threshold  value  (in  percentage  or  absolute 
values,  from  0  to  256);  all  the  pixels  with  the  brightness  exceeding  the  threshold, 
will  be  transformed  to  white,  and  those  below  the  threshold,  -  to  black  colors.  The 
choice  of  the  threshold  value  is  made  using  a  slider;  if  the  flag  of  preview  is 
herewith  set  on,  the  resulting  (scaled)  image  for  the  given  threshold  is  displayed. 

On  the  first  image  load,  the  program  evaluates  the  recommended  threshold 
value  (which  will  be  discussed  in  detail  thereafter).  According  to  our  practice,  such 
a  recommended  value  turns  out,  as  a  rule,  to  be  the  best  one. 

A  possibility  of  grayscale  to  black  and  white  image  transformation  with  the 
use  of  two  threshold  values  is  considered  in  [94].  In  doing  so,  all  shadowed  points, 
which  get  between  these  thresholds,  become  white,  while  the  others  are 
transformed  to  black.  There  is  a  number  of  alternatives  for  image  transformation, 
however,  the  choice  there  is  based  on  the  knowledge  of  specific  subject  (namely, 
which  specific  image  details  are  important;  what  is  to  be  stressed;  what  is  to  be 
removed  from  the  image,  etc.).  For  this  reason,  such  transformations  are  irrelevant 
in  the  program  aimed  at  the  generalized  determination  of  the  fractal  dimensionality 
of  images;  they  should  be  implemented  in  other  programs  for  the  operation  with 
the  graphical  information  specific  for  the  problems  of  analysis  of  images  obtained 
from  the  electron  microscopy. 
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The  button  "Evaluation  of  dimensionality"  starts  the  processing  of  the  image 
in  accordance  with  the  selected  (block  or  mass)  algorithm.  As  the  result,  one 
obtains  two  curves,  InAAv^  □  (2-Z)^)lnAr ,  or  ln(ryM(r))D  (2-D^)lnr,  as 

dependent  on  the  algorithm  used,  and  the  associated  derivative  curve  for  this 
function  (see  Fig.  83).  As  is  seen  from  the  figure,  each  plot  contains  two  curves, 
and,  in  addition,  there  is  a  switcher  "Display  the  calculation  for  the  negative".  This 
switcher  indicates  that  the  calculation  is  performed  for  two  cases.  In  the  first  case, 
the  set  P  is  represented  by  black  points  in  the  transformed  image,  and  in  the 
second  case  P  is  the  set  of  white  points.  As  a  rule,  if  the  threshold  value  is  close  to 
the  recommended  one  (which  corresponds  to  the  case  that  the  numbers  of  black 
and  white  points  are  approximately  equal),  the  curves  obtained  for  the  positive  and 
negative  turn  out  to  be  sufficiently  close  to  each  other. 

By  double-clicking  the  mouse,  the  user  can  specify  the  left  and  the  right 
boundaries  of  the  range,  where  the  numerical  value  of  the  image  fractal 
dimensionality  should  be  calculated  (i.e.,  where  the  curve  is  as  close  as  possible  to 
the  straight  line;  actually,  these  are  the  points,  where  the  derivative  reaches  its 
extremum,  and  the  second  derivative  turns  to  zero). 

The  calculated  values  are  displayed  on  the  screen,  and  the  print  button 
enables  one  to  send  the  report  to  the  printer.  In  the  last  versions  of  the  program,  at 
this  stage  the  program  estimates  the  supposed  range  for  calculations  and  offers  it 
the  user. 

While  evaluating  the  block  fractal  dimensionality,  the  question  may  arise, 
what  should  be  done  in  the  case  that  due  to  the  box  size,  the  complete  covering  of 
the  image  is  impossible  (i.e.,  a  portion  of  boxes  on  the  boundary  gets  outside  the 
image).  In  this  case,  we  employed  the  following  technique:  while  scanning  the 
image  from  the  left  to  the  right  and  from  the  top  to  the  bottom,  the  program  takes 
into  account  only  the  boxes,  which  completely  fit  into  the  image.  Next,  the 
procedure  repeats  in  the  opposite  direction,  i.e.,  from  the  right  to  the  left  and  from 
the  bottom  to  the  top.  After  that,  the  obtained  results  are  averaged. 
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The  results  obtained  in  processing  of  real  images  demonstrate,  that  such  a 
solution  for  the  above  issue  is  reasonably  justified.  In  Fig.  85,  the  curves  obtained 
for  the  same  image  of  1024x1024  pixel  size  are  shown.  On  the  left,  the  curves  for 


the  boxes  of  2',  i  =0,10  size  are  plotted,  which  completely  fit  into  the  image.  On 
the  right,  the  number  of  points  is  increased  up  to  14,  therefore,  the  boxes  are  used, 
which  get  outside  the  image.  As  is  seen  from  Fig.  85,  the  obtained  curves  are 
identical.  Therefore,  the  processing  method,  which  was  developed  and  employed  in 
the  program,  is  quite  acceptable  and  adequate  for  the  problem  under  consideration. 


Fig.  85.  Results  of  calculations  for  the  boxes  which  completely  fit  (left)  and  do  not 
fit  (right)  into  the  image. 
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Study  of  the  stability  of  the  developed  software 


As  is  known  from  practice,  the  programs  for  obtaining  the  fractal 
characteristics  do  not  always  provide  the  stable,  exact  and  reliable  results  [88].  In 
view  of  this,  we  performed  a  number  of  investigations  of  the  program  developed. 
The  need  to  have  an  a  priori  exact  information  about  the  fractal  properties  of  the 
subject  under  examination  has  pre-determined  the  Sierpinski  carpet  as  the  model 
subject  (Fig.  86)  with  the  fractal  dimensionality  Dq  -  In8/ln3  ==  1.8928 . 


Fig.86.  Sierpinski  carpet;  the  fractal  dimensionality  is  -1.893 


The  first  test  was  aimed  at  the  study  of  the  effect  of  the  number  of  included 
points  on  the  results  of  calculation.  Since  the  range  of  the  number  of  points  should 
be  in  this  case  as  wide  as  possible,  we  took  as  the  model  the  image  of  6561  x  6561 
pixel  size.  The  results  of  the  experiment  are  given  in  Fig.  87.  As  can  be  concluded 
from  the  figure,  the  effects  of  the  number  of  experimental  points  on  the 
dimensionality  are  insignificant. 


Distribution  A:  Approved  for  pubiic  reiease;  distribution  is  uniimited. 


Fig.  87.  Dependence  of  the  measured  fractal  dimensionality  on  the  number  of  the 
experimental  points. 

The  next  experiment  consisted  in  examination  of  the  effect  of  the  cyclic 
displacement  of  the  image  on  the  fractal  dimensionality.  In  this  case  we  used  the 
Sierpinski  carpet  of  smaller  size,  2187  by  2187  pixels.  As  can  be  seen  from  the 
results  of  the  numerical  experiment  given  in  Fig.  88,  the  influence  of  the 
displacement  of  the  image  on  the  calculated  fractal  dimensionality  is  rather  weak 
and  irregular. 


Fig.  88.  Dependence  of  the  measured  fractal  dimensionality  on  the  displacement  of 
the  Sierpinski  carpet. 
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For  real  images  obtained  via  electron  microscopy,  it  is  often  only  a  portion 
of  the  image  which  is  of  interest,  as,  for  instance,  in  Fig  89.  For  this  reason  one 
more  numerical  experiment  was  the  arbitrary  random  "whitening"  of  a  portion  of 
the  Sierpinski  carpet  with  subsequent  examination  of  the  effects  of  the  size  of  this 
portion  on  the  obtained  fractal  dimensionality. 


Fig.  89.  Example  with  the  subject  occupying  only  a  small  portion  of  the  image. 


In  this  case,  as  is  seen  from  Fig.  90,  we  observe  a  strong  effect  on  the 
calculated  fractal  dimensionality  for  large  areas  of  filling,  larger  than  80%. 
However,  a  closer  inspection  of  experimental  data  of  such  images  discovers  that 
there  are  two  linear  portions  on  the  curve  IniVAx^(lnAx).  As  an  example,  in  Fig. 

91,  the  experimental  data  for  Sierpinski  carpet  80%  filled  with  white  color  are 
shown.  The  fractal  dimensionality  of  Sierpinski  carpet  is  associated  with  the  curve 
of  less  steep  slope. 
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Fig.  90.  Effect  of  the  filling  a  portion  of  the  image  area  on  the  magnitude  of  the 
fractal  dimensionality. 


Fig.  91.  Experimental  data  for  the  function  In  (inAx)  for  Sierpinski  carpet 
80%  filled  with  white  color. 

Thus,  even  in  the  case  that  the  subject  to  be  examined  does  not  occupy  the 
whole  image  field,  we  can  still  determine  its  fractal  dimensionality  provided  that 
the  associated  range  of  fractality  can  be  identified. 

This  result  has  one  more  consequence  for  the  program,  which  we  elaborate. 
At  present,  we  implemented  three  versions  of  the  automatic  selection  of  the  range 
for  evaluation  of  fractality:  by  searching  the  minimum  value  for  fractal 
dimensionality;  by  removing  the  horizontal  portions  of  the  curve  InAAx^(lnAx) 
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(which  are  associated  with  fractal  dimensionality  equal  to  2);  and  by  identifying 
the  range  where  the  points  yield  the  maximum  for  the  coefficient  of  linear 
correlation.  However,  the  possible  presence  of  several  linear  portions  in  the 
experimental  function  IniVAx^  (InAv)  makes  practically  impossible  the  complete 


automation  of  the  computation  of  fractal  dimensionality  without  controlling  it  by 
user.  In  view  of  this,  the  option  for  selection  of  the  automated  determination  of  the 
fractality  range  is  provided  in  the  program. 

One  more  idea  was  inspired  by  the  operation  with  the  images  similar  to  the 
one  shown  in  Fig.  89.  It  concerns  the  possibility  of  compiling  the  image  with  the 
minimized  non-informative  component,  even  at  the  cost  of  moderately  reducing 
the  image  size.  To  simulate  such  a  treatment  with  Sierpinski  carpet,  we  performed 
the  following  numerical  experiment.  The  initial  Sierpinski  carpet  of  6561  by  6561 
pixel  size  has  been  sectioned  into  the  parts  of  a  given  size,  which  were  then  mixed 
up  at  random  by  using  the  algorithm  presented  in  [95].  After  that,  the  evaluation  of 
the  fractal  dimensionality  and  the  range  of  fractality  was  performed.  The  results  of 
this  experiment  are  given  in  Fig.  92. 


Fig.  92.  The  results  of  the  experiment  with  the  random  mixing  of  the  blocks  of  the 
Sierpinsli  carpet. 
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As  one  would  guess  a  priori,  such  a  "mixed  up"  image  should  lose  the 
information  about  its  self-similarity  with  reducing  the  block  size;  this  is  fully 
confirmed  by  the  experiment,  which  suggests  that  the  upper  limit  of  the  range  of 
image  fractality  is  with  good  accuracy  proportional  to  the  logarithm  of  block  size. 
As  for  the  fractal  dimensionality,  it  considerably  depends  on  the  block  size  as  well, 
and,  in  the  case  of  small  blocks,  is  much  less  than  its  theoretical  value.  Therefore 
the  above  treatment  based  on  putting  together  all  the  subjects  to  be  examined  in 
one  image  turns  out  to  be  a  poor  solution  to  the  problem,  since  it  can  influence  the 
subject  characteristics  to  be  evaluated. 

Further  experiments  similar  to  those  performed  in  the  study  of  the  program 
MFRDrom  [88]  are  aimed  at  the  elucidation  of  the  effects  of  size  and  orientation  of 
Sierpinski  carpet  on  its  fractal  dimensionality. 

The  effect  of  the  size  of  Sierpinski  carpet  on  its  fractal  dimensionality 
obtained  in  computations  is  demonstrated  in  Fig.  93.  It  should  be  noted  that,  along 
with  size,  the  number  of  experimental  points  are  of  importance.  For  example,  it  is 
physically  impossible  to  obtain  more  than  27  experimental  points  for  Sierpinski 
carpet  with  27  pixels  on  a  side.  Actually,  for  this  size  the  non-horizontal  portion  of 
the  curve  contains  8  points.  The  change  of  the  range  in  one  point  only  considerably 
influences  the  fractal  dimensionality  evaluated.  According  to  our  estimates,  more 
or  less  reliable  characteristics  can  be  obtained  only  for  the  images  with  at  least  500 
pixels  on  a  side. 
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Fig.  93.  Effect  of  the  size  of  Sierpinski  carpet  on  the  calculated  fractal 
dimensionality. 


Fig.  94.  Effect  of  orientation  of  Sierpinski  carpet  on  the  calculated  fractal 
dimensionality. 

Effect  of  orientation  of  Sierpinski  carpet  on  the  calculated  fractal 
dimensionality  is  demonstrated  in  Fig.  94.  As  can  be  seen,  the  rotation,  like  other 
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distortions,  which  leave  the  structure  of  Sierpinski  carpet  invariant,  influence  the 
fractal  dimensionality  insignificantly  and  in  a  irregular  manner. 

Thus,  the  numerical  experiments  performed  with  the  model  subject  have 
demonstrated  the  stable  functioning  of  the  program  and  the  stability  of  the  results, 
which  it  provides.  This  suggests  that  the  results  obtained  by  processing  the  real 
micrographs,  are  adequate. 


Developing  the  algorithm  for  obtaining  the  multifractal  characteristics  of  two- 
dimensional  images 

In  spite  of  the  success  achieved  with  employment  of  the  concept  of  regular 
fractals  for  the  evaluation  of  fractal  dimensionality  of  the  structures  typical  for 
material  science  aimed  at  the  computer  analysis,  there  is  a  number  of  studies, 
which  point  out  the  insufficiency  of  such  an  approach  [96].  Real  material  structures 
are  essentially  different  from  the  regular  fractals,  and  the  fractal  dimensionality  as 
a  numerical  characteristic  is  not  always  adequate  for  the  description  of  the 
diversity  of  the  information  contained  in  digital  images.  It  may  fail  in  qualitative 
description  of  structure  nonuniformity,  spatial  ordering  and  other  structural 
characteristics. 

It  can  be  regarded  as  the  well  established  fact  that  the  fractal  dimensionality 
as  the  only  fractal  characteristic  is  insufficient  for  the  description  of  self-similarity 
of  structures  of  real  materials.  In  the  number  of  works  (see,  for  instance,  [97,98])  it 
has  been  shown  that  such  a  possibility  can  be  provided  by  the  multifractal 
formalism. 

The  basis  of  the  multifractal  approach  for  qualitative  description  of  two- 
dimensional  structure  is  the  methodology  for  constructing  the  measure  of  set, 
which  would  approximate  the  structure  to  be  examined,  by  putting  it  in  Euclidian 
space  partitioned  into  small  square  cells  of  equal  size  (equal-cell  partition). 
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Each  cell  acquires  some  measure  represented  by  a  positive  number 
associated  with  the  portion  of  structure  occupying  this  cell.  The  ratio  of  the  number 
of  unit  structure  elements,  which  fall  into  the  cell,  to  the  total  number  of  structure 
elements,  is  accepted  as  the  sought-for  measure. 

To  define  the  measure,  one  can  use  various  physical  characteristics  of  the 
subject  under  consideration.  For  instance,  it  can  be  the  mass  of  inclusions  in  the 
matrix  being  of  interest  for  the  researcher;  in  that  case  the  total  number  of  structure 
elements  in  the  cell  is  associated  with  its  mass,  and  the  number  of  unit  elements  to 
be  counted  is  given  by  the  mass  of  inclusions. 

However,  it  should  be  pointed  out  that  the  structures  to  be  studied  in 
material  science  are  mostly  represented  by  two-dimensional  images  (e.g.,  by 
micrographs  obtained  from  the  electron  microscopy).  In  the  present-day  computer 
technology  such  images  are  represented  in  digital  form,  i.e.,  as  the  arrays  of 
discrete  equally  sized  image  elements,  the  pixels.  The  pixel  can  be  regarded  as  the 
smallest  indivisible  image  element  which  has  three  numerical  characteristics, 
namely,  the  coordinates  specifying  the  position  of  the  pixel  in  the  two-dimensional 
image,  and  its  color. 

Therefore,  for  multifractal  image  analysis,  it  is  natural  to  introduce  the 
measure  associated  with  color  properties  of  pixels  (in  the  simplest  version,  these 
are  the  two-colored  pixels  representing  the  structure  under  the  study,  along  with  its 
matrix). 

Consider  the  general  definition  of  a  multifractal.  Let  there  be  a  fractal  object 
occupying  some  limited  area  ^  of  linear  size  L  in  Euclidean  space  of  the 
dimensionality  d  (in  the  case  of  two-dimensional  image  d  =  2).  Let  it  represent,  at 
some  stage  of  its  construction,  the  set  of  □  1  points  distributed  over  this  area  in 
a  certain  way.  We  assume  that  N  ^oo , 

Let  us  partition  the  overall  area  ^  into  the  square  cells  with  the  side  eU  L 
and  the  area  .  We  will  be  interested  in  the  occupied  cells  only,  which  contain  at 
least  one  point  of  the  fractal  object.  Let  these  cells  be  numbered  as  j  =  1,N[£)  , 
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where  N^s)  is  the  total  number  of  the  occupied  cells  (which  depends  on  the  cell 
size  £). 

Denote  the  number  of  object  points  in  the  cell  with  the  index  i  as  n- . 

Then,  the  quantity 

p-{£)  =  \im  (86) 

gives  the  probability  that  a  randomly  chosen  point  belonging  to  the  constructed  set 
is  located  in  the  cell  i,  i.e.,  it  characterizes  the  relative  population  of  cells.  The 
probability  normalization  condition  yields 

The  conventional  method  of  multifractal  analysis  is  based  on  the  notion  of 
the  generalized  partition  function  Z[q,£)  where  the  power  index  can  take  any 
values  within  the  range  -oo  <q<oo\ 

=  (88) 


The  spectrum  of  the  generalized  fractal  dimensionalities  (Renyi 
dimensinalities),  which  characterize  the  given  distribution  of  points  within  the  area 
^ ,  is  defined  by  the  relation 


(89) 


where  the  scaling  exponent  r(^)is  defined,  in  turn,  as 


liq)  =  lim 


lnZ(^,£’) 


\n£ 


(90) 


For  the  case  that  Dq=D  =  const  for  all  q,  the  set  under  examination  is  the 


monofractal  described  by  one  fractal  dimensionality  D.  If  Dq  is  some  function  of  q, 
the  set  of  points  under  consideration  represents  the  multifractal. 
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By  combining  the  above  equations,  one  can  evaluate  the  Renyi 
dimensionality  as  follows: 


IrvZiq.s) 

r(,) ^^r±zz- 1 


q-l 


q-l 


-lim- 
^  - 1 


ln£ 


(91) 


The  uncertainty  of  the  0/0  type  in  resulting  from  zero  denominator  q-l 
can  be  resolved  by  using  the  evident  relation 


z:r’A’=z.T'A 


(92) 


and  performing  the  transition  to  the  limit  ^  ^  1 : 


ZN{e) 

Pi  In  Pi 

)=i 


D,  =lim 


f^O 


In^” 


(93) 


Program  implementation  of  the  algorithm  for  obtaining  the  multifractal 
characteristics  of  two-dimensional  images 

The  complex  of  programs,  which  we  developed,  computes  the  multifractal 
characteristics  of  the  structures  represented  by  the  set  of  values  of  for  different 

values  of  q.  The  calculations  are  based  on  the  above  equations,  however,  as  this 
was  in  the  case  of  program  implementation  of  the  algorithm  for  calculation  of  the 
block  fractal  dimensionality,  we  are  again  limited  to  using  the  pixel  as  the  smallest 
structure  element.  In  view  of  this,  instead  of  calculating  the  expression  (91)  in  the 
limit  of  zero  cell  size,  we  consider  the  functional  dependence 

N{e) 
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and  find  from  the  slope  angle  of  the  linear  portion,  by  approximating  the 

obtained  points  by  the  straight  line,  and  with  the  use  of  the  least  square  method. 

This  program  is  the  development  of  the  program  elaborated  earlier  and  its 
basic  features  are  the  same,  as  is  the  view  of  the  main  program  window  shown  in 
Fig.  95. 


Fig.95.  Overall  look  of  the  main  program  window. 


After  fixing  the  linear  portion  of  the  plot,  the  program  performs  the 
calculations  and  displays  the  fractal  dimensionalities  with  the  accuracy  evaluated 
in  accordance  with  the  above  relations. 

As  is  seen  from  the  screenshot  given  in  Fig.  95,  the  test  image  (Sierpinski 
carpet  of  2187x2187  pixel  size)  yields  identical  values  for  all  fractal 
dimensionalities  D  ,  which  are  consistent  with  the  theoretical  value 

D  =  ln8/ln3~l  .89279 ,  as  one  would  expect  for  this  monofractal. 
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As  was  shown  by  preliminary  investigations  similar  to  those  described 
above,  the  present  implementation  of  the  algorithm  is  rather  stable  to  the 
distortions  of  the  test  image. 

The  problem  associated  with  certain  sensitivity  of  the  employed  algorithm  to 
the  way  of  partitioning  the  image  into  the  cells  (actually,  to  the  number  of  points) 
should  be  pointed  out.  This  concerns,  in  first  place,  the  regular  fractals.  For 
instance,  in  Fig.  96,  the  experimental  points  for  Sierpinski  carpet  for  8  (top)  and  20 
(bottom)  points  of  partitioning  are  shown,  and  in  Fig.  97,  there  are  given  the  same 
results  for  the  image  of  real  structure  obtained  by  the  electron  microscopy.  As  is 
seen  from  the  given  figures,  for  the  case  of  the  regular  fractal,  the  number  of  points 
considerably  affects  the  accuracy  of  the  evaluation  of  fractal  characteristics  due  to 
the  dispersion  of  experimental  points,  whereas  in  the  case  of  real  images  this  effect 
is  much  less  pronounced. 


Fig.  96.  Experimental  curves  for  Fig.  97.  Experimental  curves  for  the 

Sierpinski  carpet  for  8  (top)  and  20  image  of  the  real  structure  for  8 

(bottom)  points.  (top)  and  20  (bottom)  points. 
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The  matter  is  that  in  the  regular  fractal  the  ideal  self-similarity  manifests 
itself  only  for  certain  size  of  trial  cells  associated  with  the  specific  structure  of  the 
fractal  itself.  For  example,  for  the  above  Sierpinski  carpet,  such  ideal  self¬ 
similarity  is  observed  in  the  case  that  the  trial  cells  are  reduced  in  size  exactly  3 
times  every  next  step.  Exactly  this  choice  of  cells  corresponds  to  the  choice  of  8 
points  for  Sierpinski  carpet  of  2187x2187  pixel  size. 

For  the  images  of  real  structures,  the  ideal  self-similarity  is  absent.  For  this 
reason,  the  experimental  data  for  different  number  of  points  used  in  calculations 
remain  close  to  each  other. 

Nevertheless,  one  of  the  options  of  further  program  development  should  be, 
probably,  the  possibility  to  analyze  the  image  size  and  calculate  the  recommended 
number  of  points  in  the  case  that  this  size  is  an  exact  power  of  some  number.  This 
would  detect  the  possible  ideal  self-similarity  and,  and  addition,  solve  the  problem 
of  boundary  effects  associated  with  the  fact  that  the  trial  cells  cannot  exactly  cover 
the  image  to  be  investigated. 


Study  of  fractal  characteristics  of  physical  objects 

The  developed  software  was  employed  in  the  Institute  for  material  science 
problems  for  practical  investigation  of  images  of  various  subjects. 

In  particular,  the  fractal  analysis  of  the  chromium  films  deposited  under  the 
argon  on  the  silicon  substrate,  as  dependent  on  the  conditions  of  production,  has 
been  performed  [99]. 

For  example,  in  Fig.  98,  there  are  shown  four  images  of  the  surface  of 
chromium  films  obtained  by  the  deposition  under  argon  at  different  ratios  of  partial 
pressures  of  oxygen  and  argon.  The  associated  dependence  of  the  fractal 
dimensionality  of  the  image  of  chromium  film  on  the  ratio  of  partial  pressures  of 
oxygen  and  argon  is  given  in  Fig  99.  As  can  be  seen,  there  is  an  expressed 


Distribution  A:  Approved  for  pubiic  reiease;  distribution  is  uniimited. 


correlation  between  the  fractal  dimensionality  of  the  image  of  the  film  and  the 
conditions  of  its  production. 


Fig.  98.  Morphology  of  the  surface  of  chromium  films  deposited  at  different  ratios 
of  partial  pressures  of  oxygen  and  argon:  0  (a),  0.88-10“^  (b),  2.99-10“^  (c), 


Fig.  99.  Fractal  dimensionality  of  the  image  of  chromium  film  as  dependent  on  the 
ratio  of  partial  pressures  of  oxygen  and  argon  during  the  deposition. 


Four  other  images  of  the  surface  of  chromium  films  obtained  by  the 
deposition  under  argon  at  different  temperatures  are  shown  in  Fig.  100. 
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Fig.  100.  Morphology  of  the  surface  of  chromium  film  deposited  at  the 
temperature  25  (a),  230  (b),  400  (c)  and  510  (d)  °C 

The  results  of  the  investigation  of  the  fractal  dimensionality  and  hardness  of 
the  films  are  given  in  Fig.  101. 


Fig.  101.  Correlation  between  the  fractal  dimensionality  and  the  hardness  of  the 
chromium  films  deposited  at  different  temperatures. 

As  we  can  see,  in  this  case  the  correlation  between  the  fractal  dimensionality 
of  the  obtained  film  and  its  physical  property,  the  hardness,  is  observed. 

Thus,  the  obtained  results  evidence  that  the  developed  software  is  applicable 
for  obtaining  the  fractal  characteristics  of  images  of  real  objects,  and  the  fractal 
characteristics  can  represent  one  more  parameter  for  the  description  of  the  objects 
under  examination. 
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lO.Conclusion 


Summarizing  the  results  obtained  we  would  like  not  only  to  highlight  and 
present  the  main  of  those  but  point  out  possible  further  directions  of  their 
development.  The  first  project  P  273  (HOARD  #068009)  made  a  solid  background 
of  computer  modeling  of  basic  components  of  technological  process  of  directed 
solidification  of  boride-boride  composites  in  systems  LaB6-MeB2(Me=Ti,  Zr,  Hf). 
The  next  project  P  510  (HOARD  #118003)  developed  a  complex  of  mathematical 
models  which  describe  this  technological  process  and  typical  properties  of 
materials-prototypes  in  the  investigated  systems  and  brought  those  to 
implementation  in  real  codes.  On  the  base  of  these  codes,  a  spectrum  of 
computative  experiments  (at  macro-,  meso-,  and  microscales)  targeted  to  develop 
the  understanding  of  peculiarities  of  eutectic  composites  was  carried  out.  The 
results  of  these  computative  experiments  are  presented  in  the  above  report  and  they 
originate  the  following  problems  for  further  work: 

I.  Macroscale 

-To  develop  the  suggested  conception  of  inner  energy  sources  up  to  a  level 
of  theory,  which  describes  the  crystallization  process  of  melts  near  eutectic  point 

-  To  perform  complex  multiscale  modeling  of  the  process  of  formation  of 
«diffusion  zone»  for  eutectic  composites  in  refractory  systems  with  a  possibility  to 
use  this  for  the  first  time  recorded  and  thermo-dynamically  grounded  formation  of 
components,  as  controlling  component  for  formation  of  required  service  properties 
of  these  composites. 

II.  Microscale 

-  to  develop  both  the  methodology  of  construction  of  interface  between 
components  of  eutectic  compositions  and  methodology  of  investigation  of 
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influence  of  atomic  parameters  of  interface  on  macro  mechanical  characteristics  of 
composites. 

-  to  develop  apparatus  of  modeling  of  self-construction  of  atomic  clusters 
which  are  components  of  the  investigated  class  of  composites  from  dimmers  to 
nanoparticles. 

III.  Mesoscale 

-  being  based  upon  close  cooperation  with  other  participants  of  the  projects 
of  «eutectic  direction»,  to  modify  programs  of  definition  of  fractal  characteristics 
of  materials  and  construction  of  cellular  automates  up  to  standard  software  of 
general  use  by  means  of  gaining  experience  while  using  those  for  different  objects 
present  in  works  of  other  participants  (initial  codes  of  computer  programs  are  in 
the  open  public  release). 
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